Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 39473) +++ trunk/Ohana/src/relastro/Makefile (revision 39474) @@ -28,4 +28,5 @@ $(SRC)/RepairWarpMeasures.$(ARCH).o \ $(SRC)/WarpImageMaps.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o \ $(SRC)/myIndex.$(ARCH).o \ $(SRC)/ConfigInit.$(ARCH).o \ @@ -112,4 +113,5 @@ $(SRC)/RepairWarpMeasures.$(ARCH).o \ $(SRC)/WarpImageMaps.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o \ $(SRC)/myIndex.$(ARCH).o \ $(SRC)/ConfigInit.$(ARCH).o \ Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 39473) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 39474) @@ -805,5 +805,5 @@ int FindWarpGroups (void); void FreeWarpGroups (void); -int GetWarpID (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y); +int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y); myIndexType *myIndexInit (); @@ -813,2 +813,6 @@ int myIndexSetEntry (myIndexType *myIndex, int value, int entry); int myIndexGetEntry (myIndexType *myIndex, int value); + +uint64_t CreatePSPSObjectID(double ra, double dec); +uint64_t CreatePSPSStackDetectionID(int sourceID, int imageID, int detID); +uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID); Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 39473) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 39474) @@ -134,4 +134,5 @@ # endif + free (image); free_astrom_table(); } Index: trunk/Ohana/src/relastro/src/RepairWarpMeasures.c =================================================================== --- trunk/Ohana/src/relastro/src/RepairWarpMeasures.c (revision 39473) +++ trunk/Ohana/src/relastro/src/RepairWarpMeasures.c (revision 39474) @@ -9,12 +9,43 @@ myAssert (!catalog->Nmeasure || catalog->measure, "programming error"); + int NfixChipID = 0, NfixStackID = 0, NfixWarpID = 0, NfixWarpImageID = 0, NmissWarp = 0, NmissStack = 0, NbadWarp = 0; + + int onePercent = catalog->Nmeasure / 100; + for (off_t j = 0; j < catalog->Nmeasure; j++) { + if (j % onePercent == 0) fprintf (stderr, "."); + Measure *measure = &catalog->measure[j]; - // we are only going to repair warp detections - if (!isGPC1warp(measure->photcode)) { - // XXX fix extID + // repair extID for non-warps: + if (isGPC1chip(measure->photcode)) { + double mjd = ohana_sec_to_mjd (measure->t); + int ccdnum = measure->photcode % 100; + uint64_t extID = CreatePSPSDetectionID (mjd, ccdnum, measure->detID); + if (extID != measure->extID) { + measure->extID = extID; + NfixChipID ++; + } continue; } + + // repair extID for non-warps: + if (isGPC1stack(measure->photcode)) { + int im = getImageByID (measure->imageID); + if (im < 0) { + if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode); + NmissStack ++; + continue; + } + uint64_t extID = CreatePSPSStackDetectionID (35, image[im].externID, measure->detID); + if (extID != measure->extID) { + measure->extID = extID; + NfixStackID ++; + } + continue; + } + + // we are only going to repair warp detections + if (!isGPC1warp(measure->photcode)) continue; // get the associated average value: @@ -23,5 +54,9 @@ // double-check for consistency - myAssert (average->objID == measure->objID, "mismatch"); + myAssert (average->objID == measure->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measure->objID); + + // warp coordinates to confirm warp + double X = measure->Xccd; + double Y = measure->Yccd; // check if this detection is far from its correct location @@ -39,17 +74,26 @@ // skip detections which are within a small distance of the expected location + // NOTE: this actually works surprisingly well near the pole double dPos = hypot(dR,dD); - if (dPos < 2.0) continue; - - // warp coordinates to confirm warp - double X = measure->Xccd; - double Y = measure->Yccd; + if (dPos < 2.0) { + int im = getImageByID (measure->imageID); + uint64_t extID = CreatePSPSStackDetectionID (34, image[im].externID, measure->detID); + if (extID != measure->extID) { + measure->extID = extID; + NfixWarpID ++; + } + continue; + } // find the corrected warp ID: - int warpID = GetWarpID (image, measure->t, measure->photcode, Rave, Dave, X, Y); - myAssert (warpID >= 0, "error finding match"); + int warpSeq = GetWarpSeq (image, measure->t, measure->photcode, Rave, Dave, X, Y); + if (warpSeq < 0) { + if (NmissWarp < 10) fprintf (stderr, "missing warp exposure: %f %f : %d %d\n", Rave, Dave, measure->t, measure->photcode); + NmissWarp ++; + continue; + } // assert on coords being TAN? - Coords *imcoords = &image[warpID].coords; + Coords *imcoords = &image[warpSeq].coords; double Rnew, Dnew; @@ -63,17 +107,26 @@ double dDnew = 3600.0*(Dnew - Dave); - // skip detections which are within a small distance of the expected location + // detections should now be repaired; fail if not dPos = hypot(dRnew,dDnew); - if (dPos > 2.0) { - // complain. abort? - fprintf (stderr, "measurement still far from average location\n"); - abort(); + if (dPos > 5.0) { + fprintf (stderr, "measurement still far from average location: %f %f vs %f %f : %d %d\n", Rave, Dave, Rnew, Dnew, measure->t, measure->photcode); + NbadWarp ++; } + measure->imageID = image[warpSeq].imageID; measure->R = Rnew; measure->D = Dnew; - // XXX fix extID + uint64_t extID = CreatePSPSStackDetectionID (34, image[warpSeq].externID, measure->detID); + if (extID != measure->extID) { + measure->extID = extID; + NfixWarpID ++; + } + NfixWarpImageID ++; } + + fprintf (stderr, "\n"); + 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); + return (TRUE); } Index: trunk/Ohana/src/relastro/src/RepairWarps.c =================================================================== --- trunk/Ohana/src/relastro/src/RepairWarps.c (revision 39473) +++ trunk/Ohana/src/relastro/src/RepairWarps.c (revision 39474) @@ -109,7 +109,8 @@ UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); - if (VERBOSE) strextend (&command, "-v"); - if (VERBOSE2) strextend (&command, "-vv"); - if (UPDATE) strextend (&command, "-update"); + if (VERBOSE) strextend (&command, "-v"); + if (VERBOSE2) strextend (&command, "-vv"); + if (UPDATE) strextend (&command, "-update"); + if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search"); fprintf (stderr, "command: %s\n", command); Index: trunk/Ohana/src/relastro/src/WarpImageMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/WarpImageMaps.c (revision 39473) +++ trunk/Ohana/src/relastro/src/WarpImageMaps.c (revision 39474) @@ -15,4 +15,5 @@ double *Dmin; double *Dmax; + int *onBoundary; } WarpGroup; @@ -120,4 +121,5 @@ ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps); ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps); + ALLOCATE (warpgroup[i].onBoundary, int, warpgroup[i].Nwarps); for (int j = 0; j < warpgroup[i].Nwarps; j++) { @@ -128,6 +130,9 @@ int Ny = image[N].NY; + warpgroup[i].onBoundary[j] = FALSE; + double R, D; XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords); + R = ohana_normalize_angle_to_midpoint(R, 180.0); warpgroup[i].Rmin[j] = R; warpgroup[i].Rmax[j] = R; @@ -137,4 +142,5 @@ // XXX need to worry about 0,360 boundary XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords); + R = ohana_normalize_angle_to_midpoint(R, 180.0); warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R); warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R); @@ -144,4 +150,5 @@ // XXX need to worry about 0,360 boundary XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords); + R = ohana_normalize_angle_to_midpoint(R, 180.0); warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R); warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R); @@ -151,9 +158,37 @@ // XXX need to worry about 0,360 boundary XY_to_RD (&R, &D, Nx, Ny, &image[N].coords); + R = ohana_normalize_angle_to_midpoint(R, 180.0); warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R); warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R); warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D); warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D); - // bump Rmin,Rmax,Dmin,Dmax by some padding? + + // bump Rmin,Rmax,Dmin,Dmax 10 arcsec worth of padding + double dR = 10.0/3600.0 / cos (RAD_DEG*warpgroup[i].Dmin[j]); + double dD = 10.0/3600.0; + + // at north pole, force test of all nearby skycells + if (warpgroup[i].Dmax[j] > 89.8) { + warpgroup[i].Rmin[j] = 0.0; + warpgroup[i].Rmax[j] = 360.0; + warpgroup[i].Dmax[j] = 90.0; + warpgroup[i].Dmin[j] -= dD; + continue; + } + + if (warpgroup[i].Rmax[j] - warpgroup[i].Rmin[j] > 270.0) { + // Rmin and Rmax are in the range 0 - 360. For images at the 0,360 boundary, + // "Rmax" is the lower edge, and "Rmin" is the upper edge. we need to flip them + // and then break the 0-360 range: + double tmp = warpgroup[i].Rmin[j]; + warpgroup[i].Rmin[j] = warpgroup[i].Rmax[j] - 360.0; + warpgroup[i].Rmax[j] = tmp; + warpgroup[i].onBoundary[j] = TRUE; + } + + warpgroup[i].Rmin[j] -= dR; + warpgroup[i].Rmax[j] += dR; + warpgroup[i].Dmin[j] -= dD; + warpgroup[i].Dmax[j] += dD; } } @@ -165,5 +200,5 @@ } -int GetWarpID (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) { +int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) { int seq = myIndexGetEntry (warpObstimeIndex, obstime); @@ -171,8 +206,17 @@ myAssert (warpgroup[seq].photcode == photcode, "oops"); + double dPosMin = NAN; + int nPosMin = -1; + // we now have the warp group, but which is the correct warp? for (int i = 0; i < warpgroup[seq].Nwarps; i++) { - if (Rave < warpgroup[seq].Rmin[i]) continue; - if (Rave > warpgroup[seq].Rmax[i]) continue; + if (warpgroup[seq].onBoundary[i]) { + int inRange1 = (Rave >= warpgroup[seq].Rmin[i]) && (Rave <= warpgroup[seq].Rmax[i]); + int inRange2 = (Rave >= warpgroup[seq].Rmin[i] + 360.0) && (Rave <= warpgroup[seq].Rmax[i] + 360.0); + if (!inRange1 && !inRange2) continue; + } else { + if (Rave < warpgroup[seq].Rmin[i]) continue; + if (Rave > warpgroup[seq].Rmax[i]) continue; + } if (Dave < warpgroup[seq].Dmin[i]) continue; if (Dave > warpgroup[seq].Dmax[i]) continue; @@ -181,20 +225,35 @@ int N = warpgroup[seq].warps[i]; - double Rwrp, Dwrp; - XY_to_RD (&Rwrp, &Dwrp, X, Y, &image[N].coords); - - float csdec = cos(Dave * RAD_DEG); + // project Rave,Dave to image pixels + double Xtst, Ytst; + RD_to_XY (&Xtst, &Ytst, Rave, Dave, &image[N].coords); + + // find the pixel offset + double dX = (Xtst - X); + double dY = (Ytst - Y); - // find the ra,dec displacement in arcsec: - double dR = 3600.0*(Rwrp - Rave)*csdec; - double dD = 3600.0*(Dwrp - Dave); - - // skip detections which are within a small distance of the expected location - double dPos = hypot(dR,dD); - if (dPos < 2.0) return N; + // skip detections which are within a small distance from the expected location + double dPos = hypot(dX,dY); + if (dPos < 20.0) return N; // 20.0 pixels = 5.0 arcsec // check for multiple possible matches?? - } - return -1; + + // record the min pos so we can see how things fail + if (isnan(dPosMin)) { + dPosMin = dPos; + nPosMin = i; + } else { + if (dPos < dPosMin) { + dPosMin = dPos; + nPosMin = i; + } + } + } + + // if we did not find any matched, give up + if (nPosMin < 0) return -1; + + // return closest image + return warpgroup[seq].warps[nPosMin]; } @@ -207,4 +266,5 @@ free (warpgroup[i].Dmin); free (warpgroup[i].Dmax); + free (warpgroup[i].onBoundary); } Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 39473) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 39474) @@ -44,5 +44,5 @@ /* unlock, if we can (else, unlocked below) and free, if we can */ - int unlockImages = !UPDATE || (RELASTRO_OP == OP_UPDATE_OFFSETS); + int unlockImages = !UPDATE || (RELASTRO_OP == OP_UPDATE_OFFSETS) || (RELASTRO_OP == OP_REPAIR_WARPS); if (unlockImages) { if (subset != image) { Index: trunk/Ohana/src/relastro/src/psps_ids.c =================================================================== --- trunk/Ohana/src/relastro/src/psps_ids.c (revision 39474) +++ trunk/Ohana/src/relastro/src/psps_ids.c (revision 39474) @@ -0,0 +1,43 @@ +# include "relastro.h" + +// Compute IDs for PSPS based on recipe they defined + +uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID) { + // t0 for detection id is 2007-01-01 00:00:00 utc + double t0 = 54101.0; + double diff = floor( 100000. * (tobs - t0) ); + int itmp = diff; + + // ccdid must be < 100 + uint64_t detectid = 1000000000*((uint64_t) itmp) + 10000000 * ((uint64_t) ccdid) + + ((uint64_t) detID); + + return detectid; +} + +uint64_t CreatePSPSStackDetectionID(int sourceID, int imageID, int detID) { + // sourceID : ID of database + table that tracked the image (< 0x7f = 127) + // imageID : external ID of the image which provided the detections (< 0x1000.0000 ~ 2.7e8) + // detID : detection sequence in image (< 0x1000.0000 ~ 2.7e8) + + assert (detID < 0x10000000); + assert (imageID < 0x10000000); + assert (sourceID < 0x7f); + + uint64_t detectid = ((uint64_t)sourceID << 56) + ((uint64_t)imageID << 28) + (uint64_t)detID; + return detectid; +} + +uint64_t CreatePSPSObjectID(double ra, double dec) { + double zh = 0.0083333; + double zid = (dec + 90.) / zh; // 0 - 180*60*2 = 21600 < 15 bits + int izone = (int) floor(zid); + double zresid = zid - ((float) izone); // 0.0 - 1.0 + + uint64_t part1, part2, part3; + part1 = (uint64_t)( izone * 10000000000000LL) ; + part2 = ((uint64_t)(ra * 1000000.)) * 10000 ; // 0 - 360*1e6 = 3.6e8 (< 29 bits) + part3 = (int) (zresid * 10000.0) ; // 0 - 10000 (1 bit == 30/10000 arcsec) (< 14 bits) + + return part1 + part2 + part3; +} Index: trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- trunk/Ohana/src/relastro/src/select_images.c (revision 39473) +++ trunk/Ohana/src/relastro/src/select_images.c (revision 39474) @@ -110,9 +110,13 @@ for (i = 0; i < Ntimage; i++) { - if (FALSE && !strncmp(timage[i].name, "o5479g0238o", 10)) { + if (FALSE && !strncmp(timage[i].name, "o6406g0242o", 10)) { fprintf (stderr, "test image 1\n"); } if (FALSE && !strncmp(timage[i].name, "o6227g0311o", 10)) { fprintf (stderr, "test image 2\n"); + } + + if (FALSE && (i >= 42819857) && (i < 42819972)) { + fprintf (stderr, "test image %s\n", timage[i].name); } @@ -173,5 +177,5 @@ DmaxImage = MAX(DmaxImage, Di[j]); } - if (RmaxImage - RminImage > 180.0) { + if (FALSE && (RmaxImage - RminImage > 180.0)) { double tmp = RminImage; RmaxImage = RminImage;