Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 39618) +++ trunk/Ohana/src/relastro/Makefile (revision 39619) @@ -25,4 +25,7 @@ RELASTRO = \ +$(SRC)/CheckMeasureToImage.$(ARCH).o \ +$(SRC)/RepairStackMeasures.$(ARCH).o \ +$(SRC)/StackImageMaps.$(ARCH).o \ $(SRC)/RepairWarps.$(ARCH).o \ $(SRC)/RepairWarpMeasures.$(ARCH).o \ @@ -110,4 +113,7 @@ RELASTRO_CLIENT = \ +$(SRC)/CheckMeasureToImage.$(ARCH).o \ +$(SRC)/RepairStackMeasures.$(ARCH).o \ +$(SRC)/StackImageMaps.$(ARCH).o \ $(SRC)/RepairWarps.$(ARCH).o \ $(SRC)/RepairWarpMeasures.$(ARCH).o \ Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 39618) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 39619) @@ -257,4 +257,16 @@ typedef struct { + int NStackFixExtID; + int NStackFixImageID; + int NstackBadCoords; + int NstackFailRepair_v1; + int NstackFailRepair_v2; + int NstackFixCoords; + int NstackMissImage; + int NstackMissTime; + int NstackNoImageID; +} StackRepairResult; + +typedef struct { int NfixChipID; int NfixStackID; @@ -269,4 +281,19 @@ } WarpRepairResult; +typedef struct { + int NstackNoImageID; + int NstackBadCoords; + int NstackBadTime; + int NstackBadImageCoords; + int NwarpNoImageID; + int NwarpBadCoords; + int NwarpBadTime; + int NwarpBadImageCoords; + int NchipNoImageID; + int NchipBadCoords; + int NchipBadTime; + int NchipBadImageCoords; +} CheckMeasureResult; + # define ID_MEAS_OBJECT_HAS_2MASS ID_MEAS_POOR_PHOTOM @@ -333,6 +360,8 @@ int USE_ICRF_POLE; +int REPAIR_STACKS; int USE_IMAGE_COORDS_FOR_REPAIR; int USE_ALL_IMAGES; +int CHECK_MEASURE_TO_IMAGE; int RESET; @@ -834,2 +863,10 @@ uint64_t CreatePSPSStackDetectionID(int sourceID, int imageID, int detID); uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID); + +void sort_by_ra (double *R, double *D, int *I, int *S, int N); +void FreeStackGroups (void); +int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y); +int MakeStackIndex (void); +int RepairStackID (StackRepairResult *result, Image *image, double Rave, double Dave, Measure *measureB, MeasureTiny *measureT); +StackRepairResult RepairStackMeasures (Catalog *catalog); +int CheckMeasureToImage (Catalog *catalog); Index: trunk/Ohana/src/relastro/src/CheckMeasureToImage.c =================================================================== --- trunk/Ohana/src/relastro/src/CheckMeasureToImage.c (revision 39619) +++ trunk/Ohana/src/relastro/src/CheckMeasureToImage.c (revision 39619) @@ -0,0 +1,271 @@ +# include "relastro.h" + +// here are the tests I am making (for chip, warp, stack): +// imageID: measure->imageID matches a valid image +// time : measure->t matches image->tzero +// coords : measure->R,D within 3 arcsec of average->R,D +// imcoords : measure->R,D within 0.1 arcsec of [r,d](image[x,y]) + +int CheckMeasureToImage (Catalog *catalog) { + + off_t Nimage; + Image *image = getimages (&Nimage, NULL); + + myAssert (!catalog->Nmeasure || catalog->measure, "programming error"); + + // XXX define StackRepairResult + CheckMeasureResult result; + result.NstackNoImageID = 0; + result.NstackBadCoords = 0; + result.NstackBadTime = 0; + result.NstackBadImageCoords = 0; + result.NwarpNoImageID = 0; + result.NwarpBadCoords = 0; + result.NwarpBadTime = 0; + result.NwarpBadImageCoords = 0; + result.NchipNoImageID = 0; + result.NchipBadCoords = 0; + result.NchipBadTime = 0; + result.NchipBadImageCoords = 0; + + int onePercent = catalog->Nmeasure / 100; + + // I want to check the measure->imageID for each photcode type: + for (off_t j = 0; j < catalog->Nmeasure; j++) { + if (j % onePercent == 0) fprintf (stderr, "."); + + Measure *measureB = &catalog->measure[j]; + + if ((measureB->averef == 0x01a) && (measureB->detID == 0x000d469a)) { + // fprintf (stderr, "example object\n"); + } + if ((measureB->averef == 0x01e) && (measureB->detID == 0x000d487f)) { + // fprintf (stderr, "example object\n"); + } + + // get the associated average value: + int averef = measureB->averef; + Average *average = &catalog->average[averef]; + + // double-check for consistency + myAssert (average->objID == measureB->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measureB->objID); + + // stack pixel coordinates + double X = measureB->Xccd; + double Y = measureB->Yccd; + + // use average position as a reference + double Rave = average->R; + double Dave = average->D; + + Rave = ohana_normalize_angle_to_midpoint(Rave, 180.0); + Dave = ohana_normalize_angle_to_midpoint(Dave, 0.0); + + // check if the new Rstk, Dstk is different from the average R,D: + float csdec = cos(Dave * RAD_DEG); + + // things to check: time, photcode, coordinates + + // do them separately to track where the failures are: + if (isGPC1stack(measureB->photcode)) { + off_t idx = getImageByID (measureB->imageID); + if (idx == -1) { + if (result.NstackNoImageID < 10) fprintf (stderr, "missing stack imageID: %f %f : %d %d %d\n", measureB->R, measureB->D, measureB->t, measureB->photcode, measureB->imageID); + result.NstackNoImageID ++; + continue; + } + + // check if the image-derived coords and the average coords agree: + Coords *imcoords = &image[idx].coords; + + double Rstk, Dstk; + XY_to_RD (&Rstk, &Dstk, X, Y, imcoords); + + Rstk = ohana_normalize_angle_to_midpoint(Rstk, 180.0); + Dstk = ohana_normalize_angle_to_midpoint(Dstk, 0.0); + + // find the ra,dec displacement in arcsec: + double dR = 3600.0*(Rstk - Rave)*csdec; + double dD = 3600.0*(Dstk - Dave); + + // 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 > 3.0) { + if (result.NstackBadCoords < 100) fprintf (stderr, "inconsistent stack avecoords: %f %f : %f %f (%d %d) 0x%08x 0x%08x\n", Rave, Dave, dR, dD, measureB->t, measureB->photcode, average->objID, average->catID); + result.NstackBadCoords ++; + } + + // check match of obstime: + if ((measureB->t < image[idx].tzero) && (measureB->t > image[idx].tzero + image[idx].NX*image[idx].trate/5000.0)) { + if (result.NstackBadTime < 10) fprintf (stderr, "inconsistent stack times: %d vs %d for %f, %f (%d)\n", measureB->t, image[idx].tzero, Rstk, Dstk, measureB->photcode); + result.NstackBadTime ++; + } + + // check if the measure coords and image-based coords match + // check if the new Rstk, Dstk is different from the average R,D: + csdec = cos(Dstk * RAD_DEG); + + // find the ra,dec displacement in arcsec: + dR = 3600.0*(Rstk - measureB->R)*csdec; + dD = 3600.0*(Dstk - measureB->D); + + dPos = hypot(dR,dD); + if (dPos > 0.1) { + if (result.NstackBadImageCoords < 10) fprintf (stderr, "inconsistent stack coords: %f %f : %f %f (%d %d)\n", Rstk, Dstk, dR, dD, measureB->t, measureB->photcode); + result.NstackBadImageCoords ++; + } + continue; + } + + // do them separately to track where the failures are: + if (isGPC1warp(measureB->photcode)) { + off_t idx = getImageByID (measureB->imageID); + if (idx == -1) { + if (result.NwarpNoImageID < 10) fprintf (stderr, "missing warp imageID: %f %f : %d %d %d\n", measureB->R, measureB->D, measureB->t, measureB->photcode, measureB->imageID); + result.NwarpNoImageID ++; + continue; + } + + // check if the image-derived coords and the average coords agree: + Coords *imcoords = &image[idx].coords; + + double Rwrp, Dwrp; + XY_to_RD (&Rwrp, &Dwrp, X, Y, imcoords); + + Rwrp = ohana_normalize_angle_to_midpoint(Rwrp, 180.0); + Dwrp = ohana_normalize_angle_to_midpoint(Dwrp, 0.0); + + // 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 + // NOTE: this actually works surprisingly well near the pole + double dPos = hypot(dR,dD); + if (dPos > 3.0) { + if (result.NwarpBadCoords < 100) fprintf (stderr, "inconsistent warp avecoords: %f %f : %f %f (%d %d) 0x%08x 0x%08x\n", Rave, Dave, dR, dD, measureB->t, measureB->photcode, average->objID, average->catID); + result.NwarpBadCoords ++; + } + + // check match of obstime: + if ((measureB->t < image[idx].tzero) && (measureB->t > image[idx].tzero + image[idx].NX*image[idx].trate/5000.0)) { + if (result.NwarpBadTime < 10) fprintf (stderr, "inconsistent warp times: %d vs %d for %f, %f (%d)", measureB->t, image[idx].tzero, Rwrp, Dwrp, measureB->photcode); + result.NwarpBadTime ++; + } + + // check if the measure coords and image-based coords match + // check if the new Rwrp, Dwrp is different from the average R,D: + csdec = cos(Dwrp * RAD_DEG); + + // find the ra,dec displacement in arcsec: + dR = 3600.0*(Rwrp - measureB->R)*csdec; + dD = 3600.0*(Dwrp - measureB->D); + + dPos = hypot(dR,dD); + if (dPos > 0.1) { + if (result.NwarpBadImageCoords < 10) fprintf (stderr, "inconsistent warp coords: %f %f : %f %f (%d %d)\n", Rwrp, Dwrp, dR, dD, measureB->t, measureB->photcode); + result.NwarpBadImageCoords ++; + } + continue; + } + + // do them separately to track where the failures are: + if (isGPC1chip(measureB->photcode)) { + off_t idx = getImageByID (measureB->imageID); + if (idx == -1) { + if (result.NchipNoImageID < 10) fprintf (stderr, "missing chip imageID: %f %f : %d %d %d\n", measureB->R, measureB->D, measureB->t, measureB->photcode, measureB->imageID); + result.NchipNoImageID ++; + continue; + } + + // check if the image-derived coords and the average coords agree: + Coords *imcoords = &image[idx].coords; + + // chip pixel coordinates + double Xchp = measureB->Xccd; + double Ychp = measureB->Yccd; + if (isfinite(measureB->Xfix) && isfinite(measureB->Yfix)) { + float dX = measureB->Xfix - measureB->Xccd; + float dY = measureB->Yfix - measureB->Yccd; + if (hypot(dX,dY) < 2.0) { + Xchp = measureB->Xfix; + Ychp = measureB->Yfix; + } + } + + double Rchp, Dchp; + XY_to_RD (&Rchp, &Dchp, Xchp, Ychp, imcoords); + + Rchp = ohana_normalize_angle_to_midpoint(Rchp, 180.0); + Dchp = ohana_normalize_angle_to_midpoint(Dchp, 0.0); + + // find the ra,dec displacement in arcsec: + double dR = 3600.0*(Rchp - Rave)*csdec; + double dD = 3600.0*(Dchp - Dave); + + // 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 > 3.0) { + if (result.NchipBadCoords < 100) fprintf (stderr, "inconsistent chip avecoords: %f %f : %f %f (%d %d) 0x%08x 0x%08x\n", Rave, Dave, dR, dD, measureB->t, measureB->photcode, average->objID, average->catID); + result.NchipBadCoords ++; + } + + // check match of obstime: + if ((measureB->t < image[idx].tzero) && (measureB->t > image[idx].tzero + image[idx].NX*image[idx].trate/5000.0)) { + if (result.NchipBadTime < 10) fprintf (stderr, "inconsistent chip times: %d vs %d for %f, %f (%d)\n", measureB->t, image[idx].tzero, Rchp, Dchp, measureB->photcode); + result.NchipBadTime ++; + } + + // check if the measure coords and image-based coords match + // check if the new Rchp, Dchp is different from the average R,D: + csdec = cos(Dchp * RAD_DEG); + + // find the ra,dec displacement in arcsec: + dR = 3600.0*(Rchp - measureB->R)*csdec; + dD = 3600.0*(Dchp - measureB->D); + + dPos = hypot(dR,dD); + if (dPos > 0.1) { + if (result.NchipBadImageCoords < 10) fprintf (stderr, "inconsistent chip coords: %f %f : %f %f (%d %d) %d 0x%08x 0x%08x | 0x%08x 0x%08x\n", Rchp, Dchp, dR, dD, measureB->t, measureB->photcode, measureB->averef, measureB->detID, measureB->imageID, average->objID, average->catID); + result.NchipBadImageCoords ++; + } + continue; + } + } + + int Nproblem = 0; + Nproblem += result.NstackNoImageID; + Nproblem += result.NstackBadCoords; + Nproblem += result.NstackBadTime; + Nproblem += result.NstackBadImageCoords; + Nproblem += result.NwarpNoImageID; + Nproblem += result.NwarpBadCoords; + Nproblem += result.NwarpBadTime; + Nproblem += result.NwarpBadImageCoords; + Nproblem += result.NchipNoImageID; + Nproblem += result.NchipBadCoords; + Nproblem += result.NchipBadTime; + Nproblem += result.NchipBadImageCoords; + + fprintf (stderr, "\n"); + fprintf (stderr, "check measure %s : %d problems : stack ImageID %d Coords %d Time %d ImCoords %d : warp ImageID %d Coords %d Time %d ImCoords %d : chip ImageID %d Coords %d Time %d ImCoords %d\n", + catalog->filename, + Nproblem, + result.NstackNoImageID, + result.NstackBadCoords, + result.NstackBadTime, + result.NstackBadImageCoords, + result.NwarpNoImageID, + result.NwarpBadCoords, + result.NwarpBadTime, + result.NwarpBadImageCoords, + result.NchipNoImageID, + result.NchipBadCoords, + result.NchipBadTime, + result.NchipBadImageCoords + ); + + return (TRUE); +} Index: trunk/Ohana/src/relastro/src/RepairStackMeasures.c =================================================================== --- trunk/Ohana/src/relastro/src/RepairStackMeasures.c (revision 39619) +++ trunk/Ohana/src/relastro/src/RepairStackMeasures.c (revision 39619) @@ -0,0 +1,205 @@ +# include "relastro.h" + +StackRepairResult RepairStackMeasures (Catalog *catalog) { + + off_t Nimage; + Image *image = getimages (&Nimage, NULL); + + myAssert (!catalog->Nmeasure || catalog->measure, "programming error"); + + // XXX define StackRepairResult + StackRepairResult result; + result.NStackFixExtID = 0; + result.NStackFixImageID = 0; + result.NstackBadCoords = 0; + result.NstackFailRepair_v1 = 0; + result.NstackFailRepair_v2 = 0; + result.NstackFixCoords = 0; + result.NstackMissImage = 0; + result.NstackMissTime = 0; + result.NstackNoImageID = 0; + + int onePercent = catalog->Nmeasure / 100; + + for (off_t j = 0; j < catalog->Nmeasure; j++) { + if (j % onePercent == 0) fprintf (stderr, "."); + + Measure *measureB = &catalog->measure[j]; + MeasureTiny *measureT = catalog->measureT ? &catalog->measureT[j] : NULL; + + // only examine stack measurements + if (!isGPC1stack(measureB->photcode)) continue; + + // get the associated average value: + int averef = measureB->averef; + Average *average = &catalog->average[averef]; + + // double-check for consistency + myAssert (average->objID == measureB->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measureB->objID); + + // stack pixel coordinates + double X = measureB->Xccd; + double Y = measureB->Yccd; + + // use average position as a reference + double Rave = average->R; + double Dave = average->D; + + Rave = ohana_normalize_angle_to_midpoint(Rave, 180.0); + Dave = ohana_normalize_angle_to_midpoint(Dave, 0.0); + + /* we have a stack detection. X,Y are correct (never modified). */ + + off_t idx = getImageByID (measureB->imageID); + if (idx == -1) { + // ** this stack detection is not correctly associated with an image. attempt to fix: + result.NstackNoImageID ++; + + if (!RepairStackID (&result, image, Rave, Dave, measureB, measureT)) { + result.NstackFailRepair_v1 ++; + } + continue; + } + + // check if the image-derived coords and the average coords agree: + Coords *imcoords = &image[idx].coords; + + double Rstk, Dstk; + XY_to_RD (&Rstk, &Dstk, X, Y, imcoords); + + Rstk = ohana_normalize_angle_to_midpoint(Rstk, 180.0); + Dstk = ohana_normalize_angle_to_midpoint(Dstk, 0.0); + + // check if the new Rstk, Dstk is different from the average R,D: + float csdec = cos(Dave * RAD_DEG); + + // find the ra,dec displacement in arcsec: + double dR = 3600.0*(Rstk - Rave)*csdec; + double dD = 3600.0*(Dstk - Dave); + + // 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 < 3.0) { + // stack image coords agree with average coords + + // check that the current extID value is good for stackdetection: + uint64_t extID = CreatePSPSStackDetectionID (35, image[idx].externID, measureB->detID); + if (extID != measureB->extID) { + measureB->extID = extID; // only in measureB, not in measureT + result.NStackFixExtID ++; + } + + // check match of obstime: + if ((measureB->t != image[idx].tzero) && (result.NstackMissTime < 10)) { + fprintf (stderr, "inconsistent stack measurement times: %d vs %d for %f, %f (%d, %ld)", measureB->t, image[idx].tzero, Rstk, Dstk, measureB->photcode, extID); + result.NstackMissTime ++; + } + + // check if the measure coords and image-based coords match + // check if the new Rstk, Dstk is different from the average R,D: + csdec = cos(Dstk * RAD_DEG); + + // find the ra,dec displacement in arcsec: + dR = 3600.0*(Rstk - measureB->R)*csdec; + dD = 3600.0*(Dstk - measureB->D); + + // 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 > 0.05) { + measureB->R = Rstk; + measureB->D = Dstk; + if (measureT) { + measureT->R = measureB->R; + measureT->D = measureB->D; + } + result.NstackFixCoords ++; + } + continue; + } + + if (!RepairStackID (&result, image, Rave, Dave, measureB, measureT)) { + result.NstackFailRepair_v2 ++; + } + } + + fprintf (stderr, "\n"); + fprintf (stderr, "repaired stacks %s : imageID %d extID %d coords %d noImageID %d : bad coords %d : failures v1 %d v2 %d image %d time %d\n", + catalog->filename, + result.NStackFixImageID, + result.NStackFixExtID, + result.NstackFixCoords, + result.NstackNoImageID, + result.NstackBadCoords, + result.NstackFailRepair_v1, + result.NstackFailRepair_v2, + result.NstackMissImage, + result.NstackMissTime); + + return (result); +} + +int RepairStackID (StackRepairResult *result, Image *image, double Rave, double Dave, Measure *measureB, MeasureTiny *measureT) { + + // get guess image based on R,D,photcode + off_t idx = GetStackSeq (image, Rave, Dave, measureB->photcode, measureB->Xccd, measureB->Yccd); + if (idx < 0) { + // still failed, this is bad + if (result->NstackMissImage < 10) fprintf (stderr, "missing stack image: %f %f : %d %d\n", Rave, Dave, measureB->imageID, measureB->photcode); + result->NstackMissImage ++; + return FALSE; + } + + // check that this image and measure now agree + if (image[idx].tzero != measureB->t) { + // still failed, this is bad + if (result->NstackMissTime < 10) fprintf (stderr, "new stack gives wrong time: %d vs %d : %f %f : %d %d\n", measureB->t, image[idx].tzero, Rave, Dave, measureB->imageID, measureB->photcode); + result->NstackMissTime ++; + return FALSE; + } + + // XXX assert on coords being TAN? + Coords *imcoords = &image[idx].coords; + + double Rnew, Dnew; + XY_to_RD (&Rnew, &Dnew, measureB->Xccd, measureB->Yccd, imcoords); + + Rnew = ohana_normalize_angle_to_midpoint(Rnew, 180.0); + Dnew = ohana_normalize_angle_to_midpoint(Dnew, 0.0); + + float csdec = cos(Dave * RAD_DEG); + + // find the ra,dec displacement in arcsec: + double dRnew = 3600.0*(Rnew - Rave)*csdec; + double dDnew = 3600.0*(Dnew - Dave); + + // detections should now be repaired; complain if not + double dPos = hypot(dRnew,dDnew); + if (dPos > 5.0) { + fprintf (stderr, "new measurement position far from average location: %f %f vs %f %f : %d %d\n", Rave, Dave, Rnew, Dnew, measureB->t, measureB->photcode); + result->NstackBadCoords ++; + } + + // fix the image ID & coords: + measureB->imageID = image[idx].imageID; + measureB->R = Rnew; + measureB->D = Dnew; + + if (measureT) { + measureT->imageID = measureB->imageID; + measureT->R = measureB->R; + measureT->D = measureB->D; + } + + // check (and fix) extID + uint64_t extID = CreatePSPSStackDetectionID (35, image[idx].externID, measureB->detID); + if (extID != measureB->extID) { + measureB->extID = extID; // only in measureB, not in measureT + result->NStackFixExtID ++; + } + + // this detection is now fixed, right? + result->NStackFixImageID ++; + return TRUE; +} Index: trunk/Ohana/src/relastro/src/StackImageMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/StackImageMaps.c (revision 39619) +++ trunk/Ohana/src/relastro/src/StackImageMaps.c (revision 39619) @@ -0,0 +1,192 @@ +# include "relastro.h" + +// generate an index of stack image ra,dec centers sorted indexes to the image IDs + +# define D_NSTACKS 1000 +# define NGROUPS 5 + +typedef struct { + unsigned short photcode; + double *Rcenter; + double *Dcenter; + int *imageID; + int *imageSeq; + int Nstacks; + int NSTACKS; +} StackGroupType; + +// one group for each photcode +StackGroupType *stackgroup = NULL; + +int MakeStackIndex (void) { + + off_t Nimage; + Image *image = getimages (&Nimage, NULL); + + ALLOCATE (stackgroup, StackGroupType, NGROUPS); + + for (int i = 0; i < NGROUPS; i++) { + stackgroup[i].photcode = 11000 + i*100; + + stackgroup[i].Nstacks = 0; + stackgroup[i].NSTACKS = D_NSTACKS; + ALLOCATE (stackgroup[i].Rcenter, double, stackgroup[i].NSTACKS); + ALLOCATE (stackgroup[i].Dcenter, double, stackgroup[i].NSTACKS); + ALLOCATE (stackgroup[i].imageID, int, stackgroup[i].NSTACKS); + ALLOCATE (stackgroup[i].imageSeq, int, stackgroup[i].NSTACKS); + } + + // first generate a list of obstime values + for (off_t i = 0; i < Nimage; i++) { + if (!isGPC1stack(image[i].photcode)) continue; + + int Ng = -1; + switch (image[i].photcode) { + case 11000: Ng = 0; break; + case 11100: Ng = 1; break; + case 11200: Ng = 2; break; + case 11300: Ng = 3; break; + case 11400: Ng = 4; break; + default: myAbort("impossible photcode"); + } + myAssert (Ng < NGROUPS, "oops"); + + int N = stackgroup[Ng].Nstacks; + XY_to_RD (&stackgroup[Ng].Rcenter[N], &stackgroup[Ng].Dcenter[N], 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords); + stackgroup[Ng].Rcenter[N] = ohana_normalize_angle (stackgroup[Ng].Rcenter[N]); + + stackgroup[Ng].imageSeq[N] = i; + stackgroup[Ng].imageID[N] = image[i].imageID; + + stackgroup[Ng].Nstacks ++; + + if (stackgroup[Ng].Nstacks >= stackgroup[Ng].NSTACKS) { + stackgroup[Ng].NSTACKS += D_NSTACKS; + REALLOCATE (stackgroup[Ng].Rcenter, double, stackgroup[Ng].NSTACKS); + REALLOCATE (stackgroup[Ng].Dcenter, double, stackgroup[Ng].NSTACKS); + REALLOCATE (stackgroup[Ng].imageID, int, stackgroup[Ng].NSTACKS); + REALLOCATE (stackgroup[Ng].imageSeq, int, stackgroup[Ng].NSTACKS); + } + } + + // sort the list: + for (int i = 0; i < NGROUPS; i++) { + sort_by_ra (stackgroup[i].Rcenter, stackgroup[i].Dcenter, stackgroup[i].imageID, stackgroup[i].imageSeq, stackgroup[i].Nstacks); + } + + return TRUE; +} + +int GetStackSeqRange (Image *image, double Rmin, double Rmax, double Dmin, double Dmax, double Rstk, double Dstk, unsigned short photcode, float X, float Y, int *nPosMin, double *dPosMin); + +// find the stack which yields the given Rstk, Dstk for the given X,Y (and photcode) +int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y) { + + if (!stackgroup) return -1; + + // we have the stack centers; find all stacks within 0.3 degrees of this point + Rstk = ohana_normalize_angle (Rstk); + double dD = 0.3; + double dR = dD / cos(RAD_DEG*Dstk); + double Rmin = Rstk - dR; + double Rmax = Rstk + dR; + double Dmin = Dstk - dD; + double Dmax = Dstk + dD; + + double dPosMin = NAN; + int nPosMin = -1; + + // if (Rmax >= 360.0) || (Rmin <= 0.0), we need to split the range into 0.0 - Ra, Rb - 360.0 + + if (Rmax >= 360.0) { + if (GetStackSeqRange (image, 0.0, Rmax - 360.0, Dmin, Dmax, Rstk, Dstk, photcode, X, Y, &nPosMin, &dPosMin)) return nPosMin; + if (GetStackSeqRange (image, Rmin, 360.0, Dmin, Dmax, Rstk, Dstk, photcode, X, Y, &nPosMin, &dPosMin)) return nPosMin; + return nPosMin; + } + + if (Rmin <= 0.0) { + if (GetStackSeqRange (image, 0.0, Rmax, Dmin, Dmax, Rstk, Dstk, photcode, X, Y, &nPosMin, &dPosMin)) return nPosMin; + if (GetStackSeqRange (image, Rmin + 360.0, 360.0, Dmin, Dmax, Rstk, Dstk, photcode, X, Y, &nPosMin, &dPosMin)) return nPosMin; + return nPosMin; + } + + // we did not find the right image, return -1 to raise an error + if (GetStackSeqRange (image, Rmin, Rmax, Dmin, Dmax, Rstk, Dstk, photcode, X, Y, &nPosMin, &dPosMin)) return nPosMin; + return nPosMin; +} + +// search for a matching stack in the given range +int GetStackSeqRange (Image *image, double Rmin, double Rmax, double Dmin, double Dmax, double Rstk, double Dstk, unsigned short photcode, float X, float Y, int *nPosMin, double *dPosMin) { + + int Ng = (photcode / 100) % 10; + int N = ohana_bisection_double (stackgroup[Ng].Rcenter, stackgroup[Ng].Nstacks, Rmin); + + for (; N < stackgroup[Ng].Nstacks; N++) { + if (stackgroup[Ng].Dcenter[N] < Dmin) continue; + if (stackgroup[Ng].Dcenter[N] > Dmax) continue; + if (stackgroup[Ng].Rcenter[N] > Rmax) break; + + int im = stackgroup[Ng].imageSeq[N]; + + // project to image pixels + double Xtst, Ytst; + RD_to_XY (&Xtst, &Ytst, Rstk, Dstk, &image[im].coords); + + // find the pixel offset + double dX = (Xtst - X); + double dY = (Ytst - Y); + + // if dPos is small, we have the right image + double dPos = hypot(dX,dY); + if (dPos < 15.0) { // 10 pixels == 2.0 arcsec + *dPosMin = dPos; + *nPosMin = im; + return TRUE; + } + + // check for multiple possible matches?? + + // record the min pos so we can see how things fail + if (isnan(*dPosMin)) { + *dPosMin = dPos; + *nPosMin = im; + } else { + if (dPos < *dPosMin) { + *dPosMin = dPos; + *nPosMin = im; + } + } + } + return FALSE; // we did not find an absolute match +} + + +void FreeStackGroups (void) { + + if (!stackgroup) return; + + for (int i = 0; i < NGROUPS; i++) { + free (stackgroup[i].Rcenter); + free (stackgroup[i].Dcenter); + free (stackgroup[i].imageID); + free (stackgroup[i].imageSeq); + } + free (stackgroup); +} + +void sort_by_ra (double *R, double *D, int *I, int *S, int N) { + +# define SWAPFUNC(A,B){ double dtmp; int itmp; \ + dtmp = R[A]; R[A] = R[B]; R[B] = dtmp; \ + dtmp = D[A]; D[A] = D[B]; D[B] = dtmp; \ + itmp = I[A]; I[A] = I[B]; I[B] = itmp; \ + itmp = S[A]; S[A] = S[B]; S[B] = itmp; \ + } +# define COMPARE(A,B)(R[A] < R[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} Index: trunk/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 39618) +++ trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 39619) @@ -11,17 +11,18 @@ image = getimages (&Nimage, NULL); - // track measurements which are far from their original position (poor at pole) - int NoffRAori = 0; - int NoffDECori = 0; + for (i = 0; i < Ncatalog; i++) { - int NoffRAchip = 0; - int NoffRAstack = 0; - int NoffRAwarp = 0; + // track measurements which are far from their original position (poor at pole) + int NoffRAori = 0; + int NoffDECori = 0; - int NoffDECchip = 0; - int NoffDECstack = 0; - int NoffDECwarp = 0; + int NoffRAchip = 0; + int NoffRAstack = 0; + int NoffRAwarp = 0; - for (i = 0; i < Ncatalog; i++) { + int NoffDECchip = 0; + int NoffDECstack = 0; + int NoffDECwarp = 0; + myAssert (!catalog[i].Nmeasure || catalog[i].measureT, "programming error"); for (j = 0; j < catalog[i].Naverage; j++) { @@ -35,4 +36,14 @@ MeasureTiny *measureT = &catalog[i].measureT[j]; Measure *measureB = catalog[i].measure ? &catalog[i].measure[j] : NULL; + + if (measureB->averef == 0x01a) { + // fprintf (stderr, "test object\n"); + } + if ((measureB->averef == 0x01a) && (measureB->detID == 0x000d469a)) { + // fprintf (stderr, "example object\n"); + } + if ((measureB->averef == 0x01e) && (measureB->detID == 0x000d487f)) { + // fprintf (stderr, "example object\n"); + } off_t im = getImageByID (measureT->imageID); @@ -111,8 +122,6 @@ } } + if (VERBOSE) fprintf (stderr, "%s : Noff ori RA %d ( %d %d %d ), Noff ori DEC %d ( %d %d %d )\n", catalog[i].filename, NoffRAori, NoffRAchip, NoffRAstack, NoffRAwarp, NoffDECori, NoffDECchip, NoffDECstack, NoffDECwarp); } - // printNcatTotal (); - - if (VERBOSE) fprintf (stderr, "Noff ori RA %d ( %d %d %d ), Noff ori DEC %d ( %d %d %d )\n", NoffRAori, NoffRAchip, NoffRAstack, NoffRAwarp, NoffDECori, NoffDECchip, NoffDECstack, NoffDECwarp); return (TRUE); Index: trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 39618) +++ trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 39619) @@ -57,4 +57,6 @@ findImages (&catalog, 1, FALSE); + if (REPAIR_STACKS) RepairStackMeasures (&catalog); + // update the detection coordinates using the new image parameters UpdateMeasures (&catalog, 1); @@ -63,4 +65,6 @@ // if Galaxy Motion Model correction is active, do not apply at this last stage UpdateObjects (&catalog, 1, -1); + + if (CHECK_MEASURE_TO_IMAGE) CheckMeasureToImage (&catalog); free_tiny_values(&catalog); @@ -190,4 +194,7 @@ if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); + if (REPAIR_STACKS) strextend (&command, "-repair-stacks"); + if (CHECK_MEASURE_TO_IMAGE) strextend (&command, "-check-measures"); + if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); if (PHOTCODE_SKIP_LIST) strextend (&command, "-photcode %s", PHOTCODE_SKIP_LIST); Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 39618) +++ trunk/Ohana/src/relastro/src/args.c (revision 39619) @@ -56,4 +56,15 @@ remove_argument (N, &argc, argv); RELASTRO_OP = OP_REPAIR_WARPS; + } + + REPAIR_STACKS = FALSE; + if ((N = get_argument (argc, argv, "-repair-stacks"))) { + remove_argument (N, &argc, argv); + REPAIR_STACKS = TRUE; + } + CHECK_MEASURE_TO_IMAGE = FALSE; + if ((N = get_argument (argc, argv, "-check-measures"))) { + remove_argument (N, &argc, argv); + CHECK_MEASURE_TO_IMAGE = TRUE; } @@ -683,4 +694,15 @@ remove_argument (N, &argc, argv); RELASTRO_OP = OP_REPAIR_WARPS; + } + + REPAIR_STACKS = FALSE; + if ((N = get_argument (argc, argv, "-repair-stacks"))) { + remove_argument (N, &argc, argv); + REPAIR_STACKS = TRUE; + } + CHECK_MEASURE_TO_IMAGE = FALSE; + if ((N = get_argument (argc, argv, "-check-measures"))) { + remove_argument (N, &argc, argv); + CHECK_MEASURE_TO_IMAGE = TRUE; } Index: trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro.c (revision 39618) +++ trunk/Ohana/src/relastro/src/relastro.c (revision 39619) @@ -46,6 +46,10 @@ } + if (REPAIR_STACKS && !PARALLEL) MakeStackIndex (); + // iterate over catalogs to make detection coordinates consistant UpdateObjectOffsets (skylist, 0, NULL); + + if (REPAIR_STACKS && !PARALLEL) FreeStackGroups (); if (!PARALLEL) { Index: trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_client.c (revision 39618) +++ trunk/Ohana/src/relastro/src/relastro_client.c (revision 39619) @@ -90,6 +90,10 @@ load_images (&db, skylist, FALSE, USE_ALL_IMAGES); + if (REPAIR_STACKS) MakeStackIndex (); + UpdateObjectOffsets (skylist, HOST_ID, HOSTDIR); + if (REPAIR_STACKS) FreeStackGroups (); + freeImages (db.ftable.buffer); gfits_db_free (&db);