Index: /trunk/magic/remove/src/Makefile.simple =================================================================== --- /trunk/magic/remove/src/Makefile.simple (revision 20519) +++ /trunk/magic/remove/src/Makefile.simple (revision 20520) @@ -5,4 +5,5 @@ streaksastrom.o \ streaksextern.o \ + warpedpixels.o \ Line.o Index: /trunk/magic/remove/src/streaks.h =================================================================== --- /trunk/magic/remove/src/streaks.h (revision 20519) +++ /trunk/magic/remove/src/streaks.h (revision 20520) @@ -64,4 +64,5 @@ pmChip *chip; // current chip pmCell *cell; // current cell + psImage *warpedPixels; psVector *tiles; float recoveryImageValue; @@ -73,4 +74,7 @@ int numCols, int numRows); +extern void computeWarpedPixels(streakFiles *sf); +extern void streaksremoveExit(psString, int); + #define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP)) #define USE_SUPPLIED_ASTROM(_stage) CHIP_LEVEL_INPUT(_stage) Index: /trunk/magic/remove/src/streaksastrom.c =================================================================== --- /trunk/magic/remove/src/streaksastrom.c (revision 20519) +++ /trunk/magic/remove/src/streaksastrom.c (revision 20520) @@ -39,5 +39,5 @@ return false; } - } else { + } else if (md) { // The metadata is the raw header // Assumes GPC1 @@ -63,4 +63,10 @@ return false; } + } else { + // no metadata regular cell + cell_x0 = 0; + cell_y0 = 0; + xParityCell = 1; + yParityCell = 1; } Index: /trunk/magic/remove/src/streaksastrom.h =================================================================== --- /trunk/magic/remove/src/streaksastrom.h (revision 20519) +++ /trunk/magic/remove/src/streaksastrom.h (revision 20520) @@ -4,4 +4,7 @@ // Structure that encapsulates Astrometry // Note: this file should have no dependence on ipp data types +// XXX: we might be able to remove this restrction. Paul S has +// used some of our datatypes in his code +// typedef struct { // opaque pointers to psModules types @@ -19,4 +22,5 @@ +#ifndef notdef // There must be some well known type lying around that we // can use for this @@ -25,4 +29,8 @@ double y; } strkPt; +#else +// TODO: remove this typedef +typedef psPlane strkPt; +#endif extern bool skyToCell(strkPt *, strkAstrom *astrom, double ra, double dec); Index: /trunk/magic/remove/src/streaksremove.c =================================================================== --- /trunk/magic/remove/src/streaksremove.c (revision 20519) +++ /trunk/magic/remove/src/streaksremove.c (revision 20520) @@ -259,5 +259,5 @@ } - // got get the astrometry + // read in the astrometry astrometry_read(sf); } @@ -312,4 +312,7 @@ } setupAstromFromFPA(sf); + if (CHIP_LEVEL_INPUT(sf->stage)) { + computeWarpedPixels(sf); + } psElemType tileType; // Type corresponding to "long" @@ -498,4 +501,11 @@ + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist")) { ; + if (CHIP_LEVEL_INPUT(stage)) { + psError(PS_ERR_UNKNOWN, true, "-skycelllist is required for raw and chip stages\n"); + return NULL; + } + } + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", NULL)) { ; if (CHIP_LEVEL_INPUT(stage)) { Index: /trunk/magic/remove/src/warpedpixels.c =================================================================== --- /trunk/magic/remove/src/warpedpixels.c (revision 20520) +++ /trunk/magic/remove/src/warpedpixels.c (revision 20520) @@ -0,0 +1,302 @@ +#include "streaks.h" +#include "pmAstrometryWCS.h" + +static void addTouchedPixels(streakFiles *sf, psString fileName); +static void nameCorners(psPlane pt[4]); +static int xLeft(psPlane *pt, int y); +static int xRight(psPlane *pt, int y); + +void +computeWarpedPixels(streakFiles *sf) +{ + pmConfig *config = sf->config; + psRegion bounds = *pmChipPixels(sf->chip); + + int width = bounds.x1 - bounds.x0; + int height = bounds.y1 - bounds.y0; + + sf->warpedPixels = psImageAlloc(width, height, PS_TYPE_U8); + psImageInit(sf->warpedPixels, 0); + + bool status; + psArray *skycells = psMetadataLookupPtr(&status, config->arguments, "SKYCELLS"); + + int n = psArrayLength(skycells); + + for (int i=0; iwarpedPixels, 0, NULL); + psFitsClose(fits); + + exit (0); +} + +static void +addTouchedPixels(streakFiles *sf, psString fileName) +{ + // read fits header + psFits *fits = psFitsOpen(fileName, "r"); + if (!fits) { + psError(PS_ERR_IO, false, "failed to open skycell file: %s", fileName); + streaksremoveExit("", PS_EXIT_DATA_ERROR); + } + psMetadata *header = psFitsReadHeader(NULL, fits); + if (!header) { + psError(PS_ERR_IO, false, "failed to read fixts header from skycell file: %s", fileName); + streaksremoveExit("", PS_EXIT_DATA_ERROR); + } + // set up astrometry for conversion from skycell image to sky + pmAstromWCS *wcs = pmAstromWCSfromHeader(header); + if (!wcs) { + psError(PS_ERR_IO, false, "failed to read astrometry from skycell file: %s", fileName); + streaksremoveExit("", PS_EXIT_DATA_ERROR); + } + int naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1"); + int naxis2 = psMetadataLookupS32(NULL, header, "NAXIS2"); + + /* now set up our wrapper to the chip astrometry to apply to the whole chip */ + sf->astrom = streakSetAstrometry(sf->astrom, sf->inAstrom->fpa, sf->chip, false, NULL, + sf->warpedPixels->numCols, sf->warpedPixels->numRows); + + + // convert corners of skycell to sky coordinates + // compute overlap with the chip using astrometry + psPlane pt[4]; + pt[0].x = 0; + pt[0].y = 0; + + pt[1].x = naxis1 - 1; + pt[1].y = 0; + + pt[2].x = naxis1 - 1; + pt[2].y = naxis2 - 1; + + pt[3].x = 0; + pt[3].y = naxis2 - 1; + + for (int i = 0; i < 4 ; i++ ) { + psSphere sky; + + // convert corner of skycell to sky coordinates + if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) { + psError(PS_ERR_IO, false, "failed to convert pt %d of %s to sky coords: %s", fileName); + streaksremoveExit("", PS_EXIT_DATA_ERROR); + } + strkPt p; + // convert to chip coordinates + if (!skyToCell(&p, sf->astrom, sky.r, sky.d)) { + psError(PS_ERR_IO, false, "failed to convert pt %d of %s to sky coords: %s", fileName); + streaksremoveExit("", PS_EXIT_DATA_ERROR); + } + pt[i].x = p.x; + pt[i].y = p.y; + } + + // put the corners in the desired order (see comments below) + + nameCorners(pt); + psString type; + if (pt[0].y < pt[2].y ) { + type = "1"; + } else { + type = "2"; + } + printf("\nSKYCELL %s Type: %s\n", fileName, type); + for (int i=0; i<4; i++) { + printf("%f %f\n", pt[i].x, pt[i].y); + } + + // Now set the touched pixels + int ymin = fmax(0, pt[1].y); + int ymax = fmin(pt[3].y, sf->warpedPixels->numRows); + printf("\nymin: %d ymax: %d\n", ymin, ymax); + for (int y = ymin ; y < ymax; y++) { + int xleft = xLeft(pt, y); + int xright = xRight(pt, y); + + if (xleft < 0) { + xleft = 0; + } + if (xleft > sf->warpedPixels->numCols) { + xleft = sf->warpedPixels->numCols - 1; + } + if (xright < 0) { + xright = 0; + } + if (xright >= sf->warpedPixels->numCols) { + xright = sf->warpedPixels->numCols - 1; + } + printf(" y: %d xleft: %d xright: %d\n", y, xleft, xright); + + psU8 *scanline = sf->warpedPixels->data.U8[y]; + + for (int x = xleft ; x <= xright; x++) { + scanline[x] = 1; + } + } +} + +// x as a function of y for the line between two points +// Note: the caller guarentees that the y's of the two points are different +static double xOfY(psPlane *pI, psPlane *pJ, int y) +{ + double m_inv = (double) (pJ->x - pI->x) / (double) (pJ->y - pI->y); + + double x = ((double) y - pI->y) * m_inv + (double) pI->x; + + return x; +} + +static int xLeft(psPlane *pt, int y) +{ + double x_d; + if (y < pt[0].y) { + // left boundary is the line from pt 0 to pt 1 + x_d = xOfY(&pt[0], &pt[1], y); + } else { + // left boundary is the line from pt 0 to pt 3 + x_d = xOfY(&pt[0], &pt[3], y); + } + return (int) x_d; +} + +static int xRight(psPlane *pt, int y) +{ + double x_d; + if (y < pt[2].y) { + // right boundary is the line from pt 1 to pt 2 + x_d = xOfY(&pt[1], &pt[2], y); + + } else { + // right boundary is the line from pt 2 to pt 3 + x_d = xOfY(&pt[2], &pt[3], y); + } + return (int) x_d; +} + +/* + * To compute the overlap of a quadrilateral with a set of + * points we transform the 4 corners to image coordinates + * and name the 4 points of the quad + pt 0 is left most (bottom most corner) + pt 1 is bottom most (right most corner) + pt 2 is right most (top most corner) + pt 3 is top most (left most corner) + +* Then we can divide the quad into 3 regions A B and C based +* on the y coordinate + +* (in the degenerate case of a rectangle aligned to the axes +* we only have 1 region) +* In each of the three regions A, B, and C the test for whether a +* point is contained in the quad on scanline y is simply + + x >= xleft(y) x < xright(y) + +* where xleft(y) and xright(y) are the bounding lines at the given +* y coordinate. + +The following diagrams show some examples. The 4 points are labeled 0 to 3 +The three regions are labeld A B and C, the boundaries between the regions +ocur at pt0.y and pt2.y + +Example 1 + + 3 + C + ---------------2 left boundary: line 0_1 y < pt0.y + B line 0_3 y >= pt0.y + 0-------------- + A right boundary: line 1_2 y < pt2.y + 1 line 2_3 y >= pt.y +************************** +Example 2 + 3 + C + 0---------------- + B + ----------------2 + A + 1 + +************************** +Example 3 + + 3 2 + + B + + 0 1 + +(region A and C are empty in the case where the points form a rectangle) +*/ + +/* + Name the corners + The following algorithm works for points that form a quadrilateral + I think it also works for the situation where 3 points are co-linear + and we have a triangle, but that isn't important for our purposes + +*/ + +static void +nameCorners(psPlane pt[4]) +{ + psPlane pt0, pt1, pt2, pt3; + + /* find pt0 the left most (bottom most) corner */ + int imin = 0; + double min = pt[0].x; + int i; + for (i=0; i<4; i++) { + double x = pt[i].x; + if ((x < min) || + ((x == min) && (pt[i].y < pt[imin].y))) { + imin = i; + min = x; + } + } + // remve pt0 from the list + pt0 = pt[imin]; + for (i=imin; i<3; i++) { + pt[i] = pt[i+1]; + } + + // find the bottom most (right most) corner from the remaining 3 + imin = 0; + min = pt[0].y; + for (i=0; i<3; i++) { + double y = pt[i].y; + if ((y < min) || + ((y == min) && (pt[i].x > pt[imin].x)) ) { + imin = i; + min = y; + } + } + // remve pt1 from the list + pt1 = pt[imin]; + for (i=imin; i<3; i++) { + pt[i] = pt[i+1]; + } + + // now find the right most (top most) of the remaining 2 points + if ((pt[0].x > pt[1].x) || + ((pt[0].x == pt[1].x) && (pt[0].y > pt[1].y))) { + pt2 = pt[0]; + pt3 = pt[1]; + } else { + pt2 = pt[1]; + pt3 = pt[0]; + } + /* write the outputs */ + pt[0] = pt0; + pt[1] = pt1; + pt[2] = pt2; + pt[3] = pt3; +}