IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39513


Ignore:
Timestamp:
Apr 2, 2016, 10:04:11 AM (10 years ago)
Author:
eugene
Message:

adding stack measure imageID repair

Location:
trunk/Ohana/src/relphot
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/Makefile

    r39481 r39513  
    3333$(SRC)/StarOps.$(ARCH).o         \
    3434$(SRC)/WarpImageMaps.$(ARCH).o   \
     35$(SRC)/StackImageMaps.$(ARCH).o  \
    3536$(SRC)/RepairWarpMeasures.$(ARCH).o      \
    3637$(SRC)/myIndex.$(ARCH).o         \
    … …  
    8889$(SRC)/StarOps.$(ARCH).o         \
    8990$(SRC)/WarpImageMaps.$(ARCH).o   \
     91$(SRC)/StackImageMaps.$(ARCH).o  \
    9092$(SRC)/RepairWarpMeasures.$(ARCH).o      \
    9193$(SRC)/myIndex.$(ARCH).o         \

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

    r39481 r39513  
    634634uint64_t CreatePSPSStackDetectionID(int sourceID, int imageID, int detID);
    635635uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID);
     636
     637int MakeStackIndex (void);
     638void FreeStackGroups (void);
     639int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y);
     640void sort_by_ra (double *R, double *D, int *I, int *S, int N);

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

    r39482 r39513  
    1818    Measure *measure = &catalog->measure[j];
    1919
     20    // get the associated average value:
     21    int averef = measure->averef;
     22    Average *average = &catalog->average[averef];
     23
     24    // double-check for consistency
     25    myAssert (average->objID == measure->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measure->objID);
     26
    2027    // repair extID for non-warps:
    2128    if (isGPC1chip(measure->photcode))  {
    … …  
    3441      int im = getImageByID (measure->imageID);
    3542      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;
     43        // if we don't have a matching image, we are probably assigned to an image on the wrong side of the sky
     44        // reconstruct the imageID from the R,D,photcode & X,Y consistency
     45        im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd);
     46        if (im < 0) {
     47          // still failed, this is bad
     48          if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
     49          NmissStack ++;
     50          continue;
     51        }
     52        // the old imageID was wrong.  replace with the new one:
     53        measure->imageID = image[im].imageID;
     54      } else {
     55        // we have a matching image, but check that it is correct by re-projecting the pixel position
     56        double Xtst, Ytst;
     57        RD_to_XY (&Xtst, &Ytst, average->R, average->D, &image[im].coords);
     58       
     59        // find the pixel offset
     60        double dX = (Xtst - measure->Xccd);
     61        double dY = (Ytst - measure->Yccd);
     62       
     63        // if dPos is small, we have the right image
     64        double dPos = hypot(dX,dY);
     65        if (dPos > 10.0) {
     66          // if dPos is large, we don't the right imageID, reconstruct the imageID from
     67          // the R,D,photcode & X,Y consistency
     68          im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd);
     69          if (im < 0) {
     70            // still failed, this is bad
     71            if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
     72            NmissStack ++;
     73            continue;
     74          }
     75        }
    3976      }
    4077      uint64_t extID = CreatePSPSStackDetectionID (35, image[im].externID, measure->detID);
    … …  
    4885    // we are only going to repair warp detections
    4986    if (!isGPC1warp(measure->photcode)) continue;
    50 
    51     // get the associated average value:
    52     int averef = measure->averef;
    53     Average *average = &catalog->average[averef];
    54 
    55     // double-check for consistency
    56     myAssert (average->objID == measure->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measure->objID);
    5787
    5888    // warp coordinates to confirm warp

  • trunk/Ohana/src/relphot/src/StackImageMaps.c

    r39512 r39513  
    44
    55# define D_NSTACKS 1000
     6# define NGROUPS 5
    67
    78typedef struct {
    … …  
    1617
    1718// one group for each photcode
    18 StackGroupType stackgroup[5];
    19 
    20 // myIndexType *warpObstimeIndex = NULL;
     19StackGroupType *stackgroup = NULL;
    2120
    2221int MakeStackIndex (void) {
    … …  
    2524  Image *image = getimages (&Nimage, NULL);
    2625
    27   for (int i = 0; i < 5; i++) {
     26  ALLOCATE (stackgroup, StackGroupType, NGROUPS);
     27
     28  for (int i = 0; i < NGROUPS; i++) {
    2829    stackgroup[i].photcode = 11000 + i*100;
    2930
    … …  
    4950      default: myAbort("impossible photcode");
    5051    }
     52    myAssert (Ng < NGROUPS, "oops");
    5153
    5254    int N = stackgroup[Ng].Nstacks;
    53     XY_to_RD (&stackgroup[Ng].Rcenter[N], &stackgroup[Ng].Dcenter[N], 0.5*image[i].NX, 0.5*image[i].NY, image[i].coords);
     55    XY_to_RD (&stackgroup[Ng].Rcenter[N], &stackgroup[Ng].Dcenter[N], 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords);
     56    stackgroup[Ng].Rcenter[N] = ohana_normalize_angle (stackgroup[Ng].Rcenter[N]);
    5457
    5558    stackgroup[Ng].imageSeq[N] = i;
    … …  
    6770  }
    6871
    69  
    70 
    71 
    7272  // sort the list:
    73   isort (obstimes, Nobstimes);
    74 
    75   // find the unique obstimes
    76   int *uniqtimes = NULL;
    77   int *uniqcount = NULL;
    78   int Nuniqtimes = 0;
    79   ALLOCATE (uniqtimes, int, Nobstimes);
    80   ALLOCATE (uniqcount, int, Nobstimes);
    81 
    82   // generate a uniq set of obstimes
    83   for (int i = 0; i < Nobstimes; Nuniqtimes ++) {
    84     uniqtimes[Nuniqtimes] = obstimes[i];
    85     int Ndup = 0;
    86     int lastValue = uniqtimes[Nuniqtimes];
    87     while ((i < Nobstimes) && (obstimes[i] == lastValue)) {
    88       i++;
    89       Ndup ++;
    90       uniqcount[Nuniqtimes] = Ndup;
    91     }
    92   }
    93   REALLOCATE (uniqtimes, int, Nuniqtimes);
    94   REALLOCATE (uniqcount, int, Nuniqtimes);
    95  
    96   // now I need to assign all warp images to a warpgroup
    97   // I have 2 options to go from warp obstime to warpgroup:
    98   // a) use a bisection lookup [may be slow]
    99   // b) use an index on obstime [require ~5*3.14*1e7*4 bytes ~ 600MB for index]
    100 
    101   // generate an index on obstime
    102   warpObstimeIndex = myIndexInit ();
    103   warpObstimeIndex->minID = obstimes[0];
    104   warpObstimeIndex->maxID = obstimes[Nobstimes-1];
    105 
    106   myIndexSetRange (warpObstimeIndex);
    107 
    108   // generate the warp groups (create the index as we go)
    109   Nwarpgroup = Nuniqtimes;
    110   ALLOCATE (warpgroup, WarpGroup, Nwarpgroup);
    111  
    112   for (int i = 0; i < Nuniqtimes; i++) {
    113     warpgroup[i].obstime = uniqtimes[i];
    114     warpgroup[i].photcode = 0; // not yet set
    115     warpgroup[i].Nwarps = 0;
    116     warpgroup[i].NWARPS = 100;
    117     ALLOCATE (warpgroup[i].warps, off_t, warpgroup[i].NWARPS);
    118     myIndexSetEntry (warpObstimeIndex, uniqtimes[i], i);
     73  for (int i = 0; i < NGROUPS; i++) {
     74    sort_by_ra (stackgroup[i].Rcenter, stackgroup[i].Dcenter, stackgroup[i].imageID, stackgroup[i].imageSeq, stackgroup[i].Nstacks);
    11975  }
    12076
    121   // assign all warps to one of the warpgroups
    122   for (off_t i = 0; i < Nimage; i++) {
    123     if (!isGPC1warp(image[i].photcode)) continue;
    124 
    125     int seq = myIndexGetEntry (warpObstimeIndex, image[i].tzero);
    126     myAssert (warpgroup[seq].obstime == image[i].tzero, "oops");
    127 
    128     if (!warpgroup[seq].photcode) {
    129       warpgroup[seq].photcode = image[i].photcode;
    130     } else {
    131       myAssert (warpgroup[seq].photcode == image[i].photcode, "oops");
    132     }
    133 
    134     int N = warpgroup[seq].Nwarps;
    135     warpgroup[seq].warps[N] = i;
    136     warpgroup[seq].Nwarps ++;
    137    
    138     CHECK_REALLOCATE (warpgroup[seq].warps, off_t, warpgroup[seq].NWARPS, warpgroup[seq].Nwarps, 100);
    139   }
    140 
    141   // we now have the warps assigned to groups.  now we need to generate the grid data to find the warp assignments
    142   for (int i = 0; i < Nwarpgroup; i++) {
    143 
    144     myAssert (warpgroup[i].Nwarps == uniqcount[i], "failure");
    145     ALLOCATE (warpgroup[i].Rmin, double, warpgroup[i].Nwarps);
    146     ALLOCATE (warpgroup[i].Rmax, double, warpgroup[i].Nwarps);
    147     ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps);
    148     ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps);
    149     ALLOCATE (warpgroup[i].onBoundary, int, warpgroup[i].Nwarps);
    150 
    151     for (int j = 0; j < warpgroup[i].Nwarps; j++) {
    152 
    153       off_t N = warpgroup[i].warps[j];
    154      
    155       int Nx = image[N].NX;
    156       int Ny = image[N].NY;
    157 
    158       warpgroup[i].onBoundary[j] = FALSE;
    159 
    160       double R, D;
    161       XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords);
    162       R = ohana_normalize_angle_to_midpoint(R, 180.0);
    163       warpgroup[i].Rmin[j] = R;
    164       warpgroup[i].Rmax[j] = R;
    165       warpgroup[i].Dmin[j] = D;
    166       warpgroup[i].Dmax[j] = D;
    167 
    168       // XXX need to worry about 0,360 boundary
    169       XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords);
    170       R = ohana_normalize_angle_to_midpoint(R, 180.0);
    171       warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
    172       warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
    173       warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
    174       warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
    175 
    176       // XXX need to worry about 0,360 boundary
    177       XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords);
    178       R = ohana_normalize_angle_to_midpoint(R, 180.0);
    179       warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
    180       warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
    181       warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
    182       warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
    183 
    184       // XXX need to worry about 0,360 boundary
    185       XY_to_RD (&R, &D, Nx, Ny, &image[N].coords);
    186       R = ohana_normalize_angle_to_midpoint(R, 180.0);
    187       warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
    188       warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
    189       warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
    190       warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
    191 
    192       // bump Rmin,Rmax,Dmin,Dmax 10 arcsec worth of padding
    193       double dR = 10.0/3600.0 / cos (RAD_DEG*warpgroup[i].Dmin[j]);
    194       double dD = 10.0/3600.0;
    195 
    196       // at north pole, force test of all nearby skycells
    197       if (warpgroup[i].Dmax[j] > 89.8) {
    198         warpgroup[i].Rmin[j] =   0.0;
    199         warpgroup[i].Rmax[j] = 360.0;
    200         warpgroup[i].Dmax[j] =  90.0;
    201         warpgroup[i].Dmin[j] -=  dD;
    202         continue;
    203       }
    204 
    205       if (warpgroup[i].Rmax[j] - warpgroup[i].Rmin[j] > 270.0) {
    206         // Rmin and Rmax are in the range 0 - 360.  For images at the 0,360 boundary,
    207         // "Rmax" is the lower edge, and "Rmin" is the upper edge.  we need to flip them
    208         // and then break the 0-360 range:
    209         double tmp = warpgroup[i].Rmin[j];
    210         warpgroup[i].Rmin[j] = warpgroup[i].Rmax[j] - 360.0;
    211         warpgroup[i].Rmax[j] = tmp;
    212         warpgroup[i].onBoundary[j] = TRUE;
    213       }
    214 
    215       warpgroup[i].Rmin[j] -=  dR;
    216       warpgroup[i].Rmax[j] +=  dR;
    217       warpgroup[i].Dmin[j] -=  dD;
    218       warpgroup[i].Dmax[j] +=  dD;
    219     }
    220   }
    221 
    222   free (uniqtimes);
    223   free (uniqcount);
    224   free (obstimes);
    22577  return TRUE;
    22678}
    22779
    228 int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) {
     80// find the stack which yields the given Rstk, Dstk for the given X,Y (and photcode)
     81int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y) {
    22982
    230   if (!warpObstimeIndex) return -1;
     83  if (!stackgroup) return -1;
    23184
    232   int seq = myIndexGetEntry (warpObstimeIndex, obstime);
    233   myAssert (seq > -1, "ooops");
    234   myAssert (warpgroup[seq].photcode == photcode, "oops");
    235  
     85  // we have the stack centers; find all stacks within 0.3 degrees of this point
     86  Rstk = ohana_normalize_angle (Rstk);
     87  double dD = 0.3;
     88  double dR = dD / cos(RAD_DEG*Dstk);
     89  double Rmin = Rstk - dR;
     90  double Rmax = Rstk + dR;
     91  double Dmin = Dstk - dD;
     92  double Dmax = Dstk + dD;
     93
     94  int Ng = (photcode / 100) % 10;
     95  int N = ohana_bisection_double (stackgroup[Ng].Rcenter, stackgroup[Ng].Nstacks, Rmin);
     96
    23697  double dPosMin = NAN;
    237   int nPosMin = -1;
     98  int    nPosMin = -1;
    23899
    239   // we now have the warp group, but which is the correct warp?
    240   for (int i = 0; i < warpgroup[seq].Nwarps; i++) {
    241     if (warpgroup[seq].onBoundary[i]) {
    242       int inRange1 = (Rave >= warpgroup[seq].Rmin[i]) && (Rave <= warpgroup[seq].Rmax[i]);
    243       int inRange2 = (Rave >= warpgroup[seq].Rmin[i] + 360.0) && (Rave <= warpgroup[seq].Rmax[i] + 360.0);
    244       if (!inRange1 && !inRange2) continue;
    245     } else {
    246       if (Rave < warpgroup[seq].Rmin[i]) continue;
    247       if (Rave > warpgroup[seq].Rmax[i]) continue;
    248     }
    249     if (Dave < warpgroup[seq].Dmin[i]) continue;
    250     if (Dave > warpgroup[seq].Dmax[i]) continue;
     100  for (; N < stackgroup[Ng].Nstacks; N++) {
     101    if (stackgroup[Ng].Dcenter[N] < Dmin) continue;
     102    if (stackgroup[Ng].Dcenter[N] > Dmax) continue;
     103    if (stackgroup[Ng].Rcenter[N] > Rmax) break;
     104   
     105    int im = stackgroup[Ng].imageSeq[N];
    251106
    252     // this is a possible image: check the coordinates:
    253     int N = warpgroup[seq].warps[i];
    254 
    255     // project Rave,Dave to image pixels
     107    // project to image pixels
    256108    double Xtst, Ytst;
    257     RD_to_XY (&Xtst, &Ytst, Rave, Dave, &image[N].coords);
    258 
     109    RD_to_XY (&Xtst, &Ytst, Rstk, Dstk, &image[im].coords);
     110   
    259111    // find the pixel offset
    260112    double dX = (Xtst - X);
    261113    double dY = (Ytst - Y);
    262114   
    263     // skip detections which are within a small distance from the expected location
     115    // if dPos is small, we have the right image
    264116    double dPos = hypot(dX,dY);
    265     if (dPos < 20.0) return N; // 20.0 pixels = 5.0 arcsec
     117    if (dPos < 10.0) return im; // 10 pixels == 2.0 arcsec
    266118
    267119    // check for multiple possible matches??
    … …  
    270122    if (isnan(dPosMin)) {
    271123      dPosMin = dPos;
    272       nPosMin = i;
     124      nPosMin = im;
    273125    } else {
    274126      if (dPos < dPosMin) {
    275127        dPosMin = dPos;
    276         nPosMin = i;
     128        nPosMin = im;
    277129      }
    278130    }
    279131  }
    280132 
    281   // if we did not find any matched, give up
    282   if (nPosMin < 0) return -1;
    283 
    284   // return closest image
    285   return warpgroup[seq].warps[nPosMin];
     133  // we did not find the right image, return what we found (-1 or closest match)
     134  return nPosMin;
    286135}
    287136
    288 void FreeWarpGroups (void) {
     137void FreeStackGroups (void) {
    289138
    290   for (int i = 0; i < Nwarpgroup; i++) {
    291     free (warpgroup[i].warps);
    292     free (warpgroup[i].Rmin);
    293     free (warpgroup[i].Rmax);
    294     free (warpgroup[i].Dmin);
    295     free (warpgroup[i].Dmax);
    296     free (warpgroup[i].onBoundary);
     139  if (!stackgroup) return;
     140
     141  for (int i = 0; i < NGROUPS; i++) {
     142    free (stackgroup[i].Rcenter);
     143    free (stackgroup[i].Dcenter);
     144    free (stackgroup[i].imageID);
     145    free (stackgroup[i].imageSeq);
    297146  }   
     147  free (stackgroup);
     148}
    298149
    299   free (warpgroup);
    300   myIndexFree (warpObstimeIndex);
     150void sort_by_ra (double *R, double *D, int *I, int *S, int N) {
     151
     152# define SWAPFUNC(A,B){ double dtmp; int itmp;  \
     153    dtmp = R[A]; R[A] = R[B]; R[B] = dtmp;      \
     154    dtmp = D[A]; D[A] = D[B]; D[B] = dtmp;      \
     155    itmp = I[A]; I[A] = I[B]; I[B] = itmp;      \
     156    itmp = S[A]; S[A] = S[B]; S[B] = itmp;      \
     157  }
     158# define COMPARE(A,B)(R[A] < R[B])
     159 
     160  OHANA_SORT (N, COMPARE, SWAPFUNC);
     161 
     162# undef SWAPFUNC
     163# undef COMPARE
     164 
    301165}

  • trunk/Ohana/src/relphot/src/relphot_objects.c

    r39511 r39513  
    4040  if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
    4141
    42   if (REPAIR_WARPS) FindWarpGroups ();
     42  if (REPAIR_WARPS) {
     43    FindWarpGroups ();
     44    MakeStackIndex ();
     45  }
    4346
    4447  // load data from each region file, only use bright stars
Note: See TracChangeset for help on using the changeset viewer.