Index: /trunk/magic/remove/src/Line.c =================================================================== --- /trunk/magic/remove/src/Line.c (revision 20280) +++ /trunk/magic/remove/src/Line.c (revision 20280) @@ -0,0 +1,347 @@ +/* Standard includes */ + +#include + +/* psLib includes */ + +#include "psMemory.h" + +/* streakremove includes */ + +#include "Line.h" + +/** Internal routine to swap double values */ + +void SwapDouble (double* first, double* second) +{ + double temp = *first; + *first = *second; + *second = temp; +} + +/** Return if a point is inside the current line segment assuming the + point and line are colinear */ + +bool InLineSegment (Line* line, strkPt* tuple) +{ + // If the current line is not vertical, check to see if the point + // is inside the two endpoints independent of direction + + if (line->begin.x != line->end.x) + { + return (tuple->x >= line->begin.x && + tuple->x <= line->end.x) || + (tuple->x >= line->end.x && + tuple->x <= line->begin.x); + } + return (tuple->y >= line->begin.y && + tuple->y <= line->end.y) || + (tuple->y >= line->end.y && + tuple->y <= line->begin.y); +} + +/** Return true if two lines intersect within their limits and in + all cases including coplanar lines, set one or two points of + intersection + + @param[in] first first line to compare for intercept + @param[in] second second line to compare for intercept + @param[out] tuple1 first intercept point if found + @param[out] tuple2 second intercept point if lines are coplanar + @param[in] inFirstSegment true if the line intercept must be within + the bounds of the first line. False if the + first line should be considered as infinite + @param[in] inSecondSegment true if the line intercept must be within + the bounds of the first line. False if the + second line should be considered as infinite + @return the number of intercept points */ + +int LineIntercept (Line* first, Line* second, strkPt* tuple1, strkPt* tuple2, + bool inFirstSegment, bool inSecondSegment) +{ + double len1, len2, w2x, w2y, t0, t1; + double dx1 = first->end.x - first->begin.x; + double dy1 = first->end.y - first->begin.y; + double dx2 = second->end.x - second->begin.x; + double dy2 = second->end.y - second->begin.y; + double wx = first->begin.x - second->begin.x; + double wy = first->begin.y - second->begin.y; + + // Calculate the perpendicular product, dt + + double dt = dx1 * dy2 - dy1 * dx2; + + // Are the lines parallel? + + if (fabs (dt) <= 1e-8) + { + // Check to see if they are overlap + + if (dx1 * wy - dy1 * wx != 0 || dx2 * wy - dy2 * wx != 0) return 0; + + // The line are coplanar, so check to see if they are degenerate + + len1 = dx1 * dx1 + dy1 * dy1; + len2 = dx2 * dx2 + dy2 * dy2; + + // First, check to see if both lines are points and if they are + // distinct tuples + + if (len1 == 0 && len2 == 0) + { + if (first->begin.x != second->begin.x || + first->begin.y != second->begin.y) return 0; + *tuple1 = first->begin; + return 1; + } + + // Check to see if the first line is a point and is inside the + // second line + + if (len1 == 0) + { + if (InLineSegment (second, &first->begin)) + *tuple1 = first->begin; + else + return 0; + } + + // Check to see if the second line is a point and is inside the + // first line + + if (len2 == 0) + { + if (InLineSegment (first, &second->begin)) + *tuple1 = second->begin; + else + return 0; + } + + // Since both lines are coplanar and have length, search for + // overlap tuples + + w2x = first->end.x - second->begin.x; + w2y = first->end.y - second->begin.y; + + t0 = (dx2 != 0) ? wx / dx2 : wy / dy2; + t1 = (dx2 != 0) ? w2x / dx2 : w2y / dy2; + + if (t0 > t1) SwapDouble (&t0, &t1); + if ((inFirstSegment || inSecondSegment) && (t0 > 1 || t1 < 0)) + { + return 0; // No overlap + } + t0 = (t0 < 0) ? 0 : t0; // Clip to min 0 + t1 = (t1 > 1) ? 1 : t1; // Clip to max 1 + + // Set the first tuple of intersection + + tuple1->x = second->begin.x + t0 * dx2; + tuple1->y = second->begin.y + t0 * dy2; + + // Check to see if the intersection is a single tuple + + if (t0 == t1) return 1; + + // Otherwise, set the second tuple of intersection as well + + tuple2->x = second->begin.x + t1 * dx2; + tuple2->y = second->begin.y + t1 * dy2; + return 2; + } + + // The segments are skew and may intersect in a tuple. + // Get the intersect parameter for the first line + + double i1 = (dx2 * wy - dy2 * wx) / dt; + if (inFirstSegment && (i1 < 0 || i1 > 1)) return 0; + + // Get the intersect parameter for the second line + + double i2 = (dx1 * wy - dy1 * wx) / dt; + if (inSecondSegment && (i2 < 0 || i2 > 1)) return 0; + + tuple1->x = first->begin.x + i1 * dx1; + tuple1->y = first->begin.y + i1 * dy1; + return 1; +} + +/** Clip the line between (0,0) and the specified upper bound + + @param[in,out] line line to be clipped within the bounds + @param[in] numCols maximum X (columns) for the line + @param[in] numRows maximum Y (rows) for the line + @return true if line overlaps the clip boundaries */ + +bool LineClip (Line *line, int numCols, int numRows) +{ + unsigned int i, found = 0; + Line boundLine, clipLine; + strkPt tuple1, tuple2, vertices[4]; + vertices[0].x = 0; vertices[0].y = 0; + vertices[1].x = numCols; vertices[1].y = 0; + vertices[2].x = numCols; vertices[2].y = numRows; + vertices[3].x = 0; vertices[3].y = numRows; + + for (i = 0; i < 4 && found < 2; ++i) + { + boundLine.begin = vertices[i]; + boundLine.end = vertices[(i + 1) % 4]; + if (LineIntercept (line, &boundLine, &tuple1, &tuple2, false, true)) + { + if (found == 0) + { + clipLine.begin = tuple1; + ++found; + } + else if (tuple1.x != clipLine.begin.x || + tuple1.y != clipLine.begin.y) + { + clipLine.end = tuple1; + ++found; + } + } + } + + // If two endpoints are found, clip the line + + if (found > 1) + { + if (clipLine.begin.x <= clipLine.end.x) + { + line->begin = clipLine.begin; + line->end = clipLine.end; + } + else + { + line->begin = clipLine.end; + line->end = clipLine.begin; + } + } + return found > 1; +} + +/** Map a line to an image for its specified width and store as a list + of pixel positions + + @param[out] pixels list of PixelPos pointers corresponding + based on the line settings + @param[in] line Line to map to pixels */ + +void PixelsFromLine (imageStreak* pixels, Line *line) +{ + PixelPos *pixel; + double slope, xOffset, yOffset, xMid, yMid, xBegin, yBegin, xEnd, yEnd, x, y; + + // Extract the endpoints + + double x1 = line->begin.x; + double y1 = line->begin.y; + double x2 = line->end.x; + double y2 = line->end.y; + + // Determine the width and height of the line + + double dx = x2 - x1; + double dy = y2 - y1; + double dr = sqrt (dx * dx + dy * dy); + double halfWidth = line->width / 2.0; + double halfWidth2 = halfWidth * halfWidth; + if (!dr) return; + + // Step point by point based on the dominate axis + + if (fabs (dx) > fabs (dy)) + { + // --------------------- + // horizontal(ish) lines + // --------------------- + // If line is back to front, reorder endpoints and recompute + // the line width and height + + if (x1 > x2) + { + SwapDouble (&x1, &x2); + SwapDouble (&y1, &y2); + dx = x2 - x1; + dy = y2 - y1; + } + + // Compute the slope of the line + + slope = dy / dx; + + // Compute the x and y offsets for the line width extent + + xOffset = halfWidth * dy / dr; + yOffset = halfWidth * dr / dx; + yMid = y1 + slope * (floor (x1 - xOffset) - x1); + xBegin = floor (x1 - xOffset); + xEnd = ceil (x2 + xOffset) + 1.0; + + for (x = xBegin; x != xEnd; ++x) + { + yBegin = floor (yMid - yOffset); + yEnd = ceil (yMid + yOffset) + 1.0; + for (y = yBegin; y != yEnd; ++y) + { + if ((x * x + y * y) <= halfWidth2) + { + pixel = psAlloc (sizeof(PixelPos)); + pixel->x = (int) x; + pixel->y = (int) y; + psArrayAdd (pixels, 1024, pixel); + psFree (pixel); + } + } + yMid += slope; + } + } + else + { + // ------------------- + // vertical(ish) lines + // ------------------- + // If line is back to front, reorder endpoints and recompute + // the line width and height + + if (y1 > y2) + { + SwapDouble (&x1, &x2); + SwapDouble (&y1, &y2); + dx = x2 - x1; + dy = y2 - y1; + } + + // Compute the slope of the line + + slope = dx / dy; + + // Compute the x and y offsets for the line width extent + + xOffset = halfWidth * dr / dy; + yOffset = halfWidth * dx / dr; + + xMid = x1 + slope * (floor (y1 - yOffset) - y1); + yBegin = floor (y1 - yOffset); + yEnd = ceil (y2 + yOffset) + 1.0; + + for (y = yBegin; y != yEnd; ++y) + { + xBegin = floor (xMid - xOffset); + xEnd = ceil (xMid + xOffset) + 1.0; + for (x = xBegin; x != xEnd; ++x) + { + if ((x * x + y * y) <= halfWidth2) + { + pixel = psAlloc (sizeof(PixelPos)); + pixel->x = (int) x; + pixel->y = (int) y; + psArrayAdd (pixels, 1024, pixel); + psFree (pixel); + } + } + xMid += slope; + } + } // if (fabs (dx) > fabs (dy)) +} Index: /trunk/magic/remove/src/Line.h =================================================================== --- /trunk/magic/remove/src/Line.h (revision 20280) +++ /trunk/magic/remove/src/Line.h (revision 20280) @@ -0,0 +1,32 @@ +#ifndef STREAK_LINE_H +#define STREAK_LINE_H + +#include "streaksextern.h" + +/** Define a line object converted from streak equatorial object */ + +typedef struct { + strkPt begin; + strkPt end; + double width; +} Line; + +/** Clip the line between (0,0) and the specified upper bound + + @param[in,out] line line to be clipped within the bounds + @param[in] numCols maximum X (columns) for the line + @param[in] numRows maximum Y (rows) for the line + @return true if line overlaps the clip boundaries */ + +extern bool LineClip (Line *line, int numCols, int numRows); + +/** Map a line to an image for its specified width and append as + a list of pixel positions + + @param[out] pixels list of PixelPos pointers corresponding + based on the line settings + @param[in] line Line to map to pixels */ + +extern void PixelsFromLine (imageStreak* pixels, Line *line); + +#endif /* STREAK_LINE_H */ Index: /trunk/magic/remove/src/Makefile.simple =================================================================== --- /trunk/magic/remove/src/Makefile.simple (revision 20280) +++ /trunk/magic/remove/src/Makefile.simple (revision 20280) @@ -0,0 +1,27 @@ +# skeleton Makefile for streaksremove + +OBJECTS= \ + streaksremove.o \ + streaksastrom.o \ + streaksextern.o \ + Line.o + +CFLAGS=`psmodules-config --cflags` +LDFLAGS=`psmodules-config --libs` + +streaksremove: ${OBJECTS} + +clean: + rm -f *.o streaksremove + + +clean-local: + -rm -f TAGS + +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` + +### Error codes. +BUILT_SOURCES = +CLEANFILES = Index: /trunk/magic/remove/src/streaks.h =================================================================== --- /trunk/magic/remove/src/streaks.h (revision 20280) +++ /trunk/magic/remove/src/streaks.h (revision 20280) @@ -0,0 +1,83 @@ +#ifndef STREAKS_H +#define STREAKS_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "string.h" +#include "pslib.h" +#include "psmodules.h" + +typedef struct { + int dummy; +} Warps; + +#include "streaksastrom.h" + +// WES file names are ipp uris usually neb:// +// we need to resolve these to unix file names +// sFile is a wrapper containing the information for our files + +typedef struct { + psString resolved_name; + psString name; + bool inNebulous; + psFits *fits; + int nHDU; + psString extname; + psMetadata *header; + int numZPlanes; + psImage *image; + psArray *imagecube; + pmFPAfile *pmfile; + int numCols; + int numRows; +} sFile; + + +typedef enum { + IPP_STAGE_NONE = 0, + IPP_STAGE_RAW, + IPP_STAGE_CHIP, + IPP_STAGE_WARP, + IPP_STAGE_DIFF +} ippStage; + + +typedef struct { + pmConfig *config; + ippStage stage; + int extnum; + int nHDU; // number of HDUs in the inputs (all must be equal) + sFile *inImage; + sFile *inMask; + sFile *inWeight; + sFile *outImage; + sFile *outMask; + sFile *outWeight; + sFile *recImage; + sFile *recMask; + sFile *recWeight; + psString class_id; + pmFPAfile *inAstrom; + strkAstrom *astrom; + bool bilevelAstrometry; + pmFPAview *view; + pmChip *chip; // current chip + pmCell *cell; // current cell + psVector *tiles; + float recoveryImageValue; + float recoveryMaskValue; + float recoveryWeightValue; +} streakFiles; + +extern strkAstrom *streakSetAstrometry(strkAstrom *, pmFPA *, pmChip *, bool fromCell, psMetadata *md, + int numCols, int numRows); + +#define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP)) +#define USE_SUPPLIED_ASTROM(_stage) CHIP_LEVEL_INPUT(_stage) + +#define SFILE_IS_IMAGE(_sfile) (_sfile->image || _sfile->imagecube) + +#endif // STREAKS_H Index: /trunk/magic/remove/src/streaksastrom.c =================================================================== --- /trunk/magic/remove/src/streaksastrom.c (revision 20280) +++ /trunk/magic/remove/src/streaksastrom.c (revision 20280) @@ -0,0 +1,132 @@ +#include "streaks.h" + +// Initialize the astrometry object for current cell +strkAstrom * +streakSetAstrometry(strkAstrom *astrom, pmFPA *fpa, pmChip *chip, bool fromCell, psMetadata *md, int numCols, int numRows) +{ + if (!astrom) { + astrom = psAlloc(sizeof(strkAstrom)); + astrom->pt = pmAstromObjAlloc(); + } + astrom->fpa = fpa; + astrom->chip = chip; + + int cell_x0, cell_y0; + int xParityCell, yParityCell; + bool mdok; + + if (fromCell) { + // The metadata is the cell concepts + cell_x0 = psMetadataLookupS32(&mdok, md, "CELL.X0"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CELL.X0 for cell is not set.\n"); + return NULL; + } + cell_y0 = psMetadataLookupS32(&mdok, md, "CELL.Y0"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "CELL.Y0 for cell is not set.\n"); + return NULL; + } + + xParityCell = psMetadataLookupS32(&mdok, md, "CELL.XPARITY"); + if (!mdok || (xParityCell != 1 && xParityCell != -1)) { + psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY for cell is not set.\n"); + return false; + } + yParityCell = psMetadataLookupS32(&mdok, md, "CELL.YPARITY"); + if (!mdok || (yParityCell != 1 && yParityCell != -1)) { + psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY for cell is not set.\n"); + return false; + } + } else { + // The metadata is the raw header + // Assumes GPC1 + cell_x0 = psMetadataLookupS32(&mdok, md, "IMNPIX1"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "IMNPIX1 for cell is not set.\n"); + return NULL; + } + cell_y0 = psMetadataLookupS32(&mdok, md, "IMNPIX2"); + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "IMNPIX2 for cell is not set.\n"); + return NULL; + } + + xParityCell = psMetadataLookupS32(&mdok, md, "ATM1_1"); + if (!mdok || (xParityCell != 1 && xParityCell != -1)) { + psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY for cell is not set.\n"); + return false; + } + yParityCell = psMetadataLookupS32(&mdok, md, "ATM2_2"); + if (!mdok || (yParityCell != 1 && yParityCell != -1)) { + psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY for cell is not set.\n"); + return false; + } + } + + astrom->cell_x0 = (double) cell_x0; + astrom->cell_y0 = (double) cell_y0; + // XXX: make sure we're applying cell parity correctly + // the mosaic code uses (xParityCell * xParityChip == xParityTarget) ? -1 : 1 + // (yParityCell * yParityChip == yParityTarget) ? -1 : 1 + // since the transform already accounts for the chip parity we just use the cell parity + astrom->xParity = (double) xParityCell; + astrom->yParity = (double) yParityCell; + + astrom->numCols = numCols; + astrom->numRows = numRows; + + return astrom; +} + +bool +skyToCell(strkPt *outPt, strkAstrom *astrom, double ra, double dec) +{ + pmFPA *fpa = (pmFPA *) astrom->fpa; + pmChip *chip = (pmChip *) astrom->chip; + pmAstromObj *pt = (pmAstromObj *) astrom->pt; + + pt->sky->r = ra; + pt->sky->d = dec; + pt->sky->rErr = 0; // Is setting these errors to zero the right thing to do? + pt->sky->dErr = 0; + + // convert from sky to TP to FP + psProject(pt->TP, pt->sky, fpa->toSky); + psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP); + // convert from FP to chip + psPlaneTransformApply(pt->chip, chip->fromFPA, pt->FP); + + // convert from chip to cell + outPt->x = (pt->chip->x - astrom->cell_x0) * astrom->xParity; + outPt->y = (pt->chip->y - astrom->cell_y0) * astrom->yParity; + + return true; +} + +bool +cellToSky(strkPt *outPt, strkAstrom *astrom, double x, double y) +{ + pmFPA *fpa = (pmFPA *) astrom->fpa; + pmChip *chip = (pmChip *) astrom->chip; + pmAstromObj *pt = (pmAstromObj *) astrom->pt; + + // XXX: TODO: confirm that this cell to chip transformation is correct + pt->chip->x = (x * astrom->xParity) + astrom->cell_x0; + pt->chip->y = (y * astrom->yParity) + astrom->cell_y0; + + pt->chip->xErr = 0; // Is setting these errors to zero the right thing to do? + pt->chip->yErr = 0; + + // chip to FP + psPlaneTransformApply(pt->FP, chip->toFPA, pt->chip); + // FP to TP to sky + psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP); + psDeproject(pt->sky, pt->TP, fpa->toSky); + + outPt->x = pt->sky->r; + outPt->y = pt->sky->d; + + return true; +} + Index: /trunk/magic/remove/src/streaksastrom.h =================================================================== --- /trunk/magic/remove/src/streaksastrom.h (revision 20280) +++ /trunk/magic/remove/src/streaksastrom.h (revision 20280) @@ -0,0 +1,31 @@ +#ifndef STREAKSASTROM_H +#define STREAKSASTROM_H + +// Structure that encapsulates Astrometry +// Note: this file should have no dependence on ipp data types +typedef struct { + // opaque pointers to psModules types + void *pt; // pmAstromObj * + void *fpa; // pmFPA * + void *chip; // pmChip * + + double cell_x0; + double cell_y0; + double xParity; // 1 or -1 + double yParity; + int numCols; + int numRows; +} strkAstrom; + + +// There must be some well known type lying around that we +// can use for this +typedef struct { + double x; + double y; +} strkPt; + +extern bool skyToCell(strkPt *, strkAstrom *astrom, double ra, double dec); +extern bool cellToSky(strkPt *, strkAstrom *astrom, double x, double y); + +#endif // STREAKS_ASTROM_H Index: /trunk/magic/remove/src/streaksremove.c =================================================================== --- /trunk/magic/remove/src/streaksremove.c (revision 20280) +++ /trunk/magic/remove/src/streaksremove.c (revision 20280) @@ -0,0 +1,977 @@ +#include "streaks.h" +#include "streaksextern.h" + +extern bool sFileLock(sFile * sfile); +extern bool sFileUnlock(sFile * sfile); +extern bool sFileReplace(sFile *dest, sFile *src); + +static void bail_out(psString str, int exitCode) { + psErrorStackPrint(stderr, str); + exit(exitCode); +} + +static void +sFileFree(sFile *sfile) +{ + psFree(sfile->resolved_name); + psFree(sfile->name); + psFree(sfile); +} + +psString +resolveName(pmConfig *config, sFile *sfile) +{ + // TODO: do the nebulous lookup + sfile->inNebulous = false; + + return psStringCopy(sfile->name); +} + +sFile *sFileOpen(pmConfig *config, ippStage stage, psString fileSelect, + psString outputExt, bool required) +{ + bool status; + sFile *sfile = psAlloc(sizeof(sFile)); + memset(sfile, 0, sizeof(sFile)); + + // if stage is raw or chip use psFits for input, if warp or diff start with FPA + // so we have it for astrometry + if ((CHIP_LEVEL_INPUT(stage)) || strcmp(fileSelect, "INPUT")) { + sfile->name = psMetadataLookupStr(&status, config->arguments, fileSelect); + if (!status || !sfile->name) { + if (required) { + psError(PS_ERR_IO, false, "Failed to lookup name for %s", fileSelect); + sFileFree(sfile); + bail_out("", 1); + } + return NULL; + } + } else { + // stage is warp or diff AND fileSelect eq "INPUT" + // get data from pmFPAfile + sfile->pmfile = pmFPAfileDefineFromArgs(&status, config, "PPSUB.INPUT", "INPUT"); + if (!sfile->pmfile) { + psError(PS_ERR_IO, false, "Failed to defile file for name for %s", fileSelect); + bail_out("", PS_EXIT_DATA_ERROR); + } + pmFPAview *view = pmFPAviewAlloc(0); + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "Failed to load input."); + bail_out("", PS_EXIT_DATA_ERROR); + } + psFree(view); + sfile->name = sfile->pmfile->filename; + sfile->resolved_name = psStringCopy(sfile->name); + // copy header from fpu + sfile->header = (psMetadata*) psMemIncrRefCounter(sfile->pmfile->fpa->hdu->header); + return sfile; + } + + // if outputExt is not null it contains the extension to append to name + // and the file is to be opened for writing + if (outputExt) { + psStringAppend(&sfile->name, outputExt); + sfile->resolved_name = psStringCopy(sfile->name); + sfile->fits = psFitsOpen(sfile->resolved_name, "w"); + sfile->fits->options = psFitsOptionsAlloc(); + } else { + sfile->resolved_name = resolveName(config, sfile); + if (!sfile->resolved_name) { + psError(PS_ERR_IO, false, "Failed to resolve name for %s", sfile->name); + sFileFree(sfile); + bail_out("", 1); + } + sfile->fits = psFitsOpen(sfile->resolved_name, "r"); + if (sfile->fits) { + sfile->nHDU = psFitsGetSize(sfile->fits); + } + } + + if (!sfile->fits) { + psError(PS_ERR_IO, false, "failed to open fits file %s for %s", + sfile->resolved_name, outputExt ? "writing" : "reading"); + sFileFree(sfile); + bail_out("", 1); + } + + return sfile; +} + +#ifdef notyet + sfile->file = pmFPAfileDefineFromArgs(&status, config, name, fileSelect); + if (!status || !sfile->file) { + if (required) { + psError(PS_ERR_IO, false, "Failed to build FPA for %s", name); + } + return NULL; + } +#endif + + +bool +astrometry_read(streakFiles *sf) +{ + pmHDU *phu = pmFPAviewThisPHU(sf->view, sf->inAstrom->fpa); + if (phu) { + char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + if (ctype) { + sf->bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } + } + + if (sf->bilevelAstrometry) { +#ifdef notyet + // Do we get here for GPC1 ? I don't necessarily have an fpa for the input image + if (!pmAstromReadBilevelMosaic(sf->inImage->file->fpa, phu->header)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA."); + return false; + } +#endif + } else { + pmHDU *hdu = pmFPAviewThisHDU(sf->view, sf->inAstrom->fpa); + PS_ASSERT_PTR_NON_NULL(hdu, 1) + PS_ASSERT_PTR_NON_NULL(hdu->header, 1) + + // we use a default FPA pixel scale of 1.0 + if (!pmAstromReadWCS (sf->inAstrom->fpa, sf->chip, hdu->header, 1.0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA."); + return false; + } + } + + return true; +} + +static void +setupAstromFromFPA(streakFiles *sf) +{ + sf->view = pmFPAviewAlloc(0); + + if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "Failed to load input."); + bail_out("", PS_EXIT_DATA_ERROR); + } + + while ((sf->chip = pmFPAviewNextChip(sf->view, sf->inAstrom->fpa, 1))) { + if (sf->inAstrom->fpa->chips->n == 1) { + // There's only one chip in this FPA and we've got it so break + break; + } + bool status; + psString chip_name = psMetadataLookupStr(&status, sf->chip->concepts, "CHIP.NAME"); + if (!strcmp(chip_name, sf->class_id)) { + break; + } + } + if (!sf->chip) { + psError(PS_ERR_UNKNOWN, true, "Failed to find chip with data."); + bail_out("", PS_EXIT_DATA_ERROR); + } + if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed to load chip"); + bail_out("", PS_EXIT_DATA_ERROR); + } + + // got get the astrometry + astrometry_read(sf); +} + +streakFiles *openFiles(pmConfig *config) +{ + bool status; + streakFiles *sf = psAlloc(sizeof(streakFiles)); + memset(sf, 0, sizeof(*sf)); + + sf->config = config; + + // error checking is done by sFileOpen. If a file can't be open we + // exit + ippStage stage = psMetadataLookupS32(&status, config->arguments, "STAGE"); + + sf->stage = stage; + sf->extnum = 0; + + sf->class_id = psMetadataLookupStr(&status, config->arguments, "CLASS_ID"); + + sf->inImage = sFileOpen(config, stage, "INPUT", NULL, true); + sf->nHDU = sf->inImage->nHDU; + sf->outImage = sFileOpen(config, stage, "OUTPUT", ".fits", true); + sf->recImage = sFileOpen(config, stage, "RECOVERY", ".fits", true); + + sf->inMask = sFileOpen(config, stage, "INPUT.MASK", NULL, false); + if (sf->inMask) { + sf->outMask = sFileOpen(config, stage, "OUTPUT", ".mk.fits", true); + sf->recMask = sFileOpen(config, stage, "RECOVERY", "wt.fits", true); + } + sf->inWeight = sFileOpen(config, stage, "INPUT.WEIGHT", NULL, false); + if (sf->inWeight) { + sf->outWeight = sFileOpen(config, stage, "OUTPUT", ".wt.fits", true); + sf->recWeight = sFileOpen(config, stage, "RECOVERY", ".wt.fits", true); + } + + // load astrometry if supplied + if (USE_SUPPLIED_ASTROM(sf->stage)) { + sf->inAstrom = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM", "ASTROM"); + } else { + sf->inAstrom = sf->inImage->pmfile; + if (!sf->inImage->pmfile) { + bail_out("unexpected null pmFPAfile", PS_EXIT_CONFIG_ERROR); + } + } + setupAstromFromFPA(sf); + + psElemType tileType; // Type corresponding to "long" + if (sizeof(long) == sizeof(psS64)) { + tileType = PS_TYPE_S64; + } else if (sizeof(long) == sizeof(psS32)) { + tileType = PS_TYPE_S32; + } else { + psAbort("can't map (long) type to a psLib type"); + } + + sf->tiles = psVectorAlloc(3, tileType); // Tile sizes + if (tileType == PS_TYPE_S64) { + sf->tiles->data.S64[0] = 0; + sf->tiles->data.S64[1] = 1; + sf->tiles->data.S64[2] = 1; + } else { + sf->tiles->data.S32[0] = 0; + sf->tiles->data.S32[1] = 1; + sf->tiles->data.S32[2] = 1; + } + + return sf; +} + + +bool +streakFilesLock(streakFiles *sf) +{ + if (sf->inImage->inNebulous ) { + // TODO: make the nebulous call to lock the file + } + return true; +} +bool +streakFilesUnlock(streakFiles *sf) +{ + if (sf->inImage->inNebulous ) { + // TODO: make the nebulous call to unlock the file + } + return true; +} + +static bool +moveExt(sFile *sfile, int extnum) +{ + bool status = psFitsMoveExtNum(sfile->fits, extnum, false); + if (!status) { + psError(PS_ERR_IO, false, + "failed to move to extension %d for %s", extnum, sfile->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + return true; +} + +static bool +streakFilesNextExtension(streakFiles *sf) +{ + // unless stage is raw or chip when we get here we're done + if (!CHIP_LEVEL_INPUT(sf->stage)) { + return false; + } + + ++sf->extnum; + if (sf->nHDU == 1) { + // return true the first time through + return (sf->extnum == 1); + } + if (sf->extnum < sf->nHDU) { + moveExt(sf->inImage, sf->extnum); + if (sf->inMask) { + moveExt(sf->inMask, sf->extnum); + } + if (sf->inWeight) { + moveExt(sf->inWeight, sf->extnum); + } + return true; + } else { + return false; + } +} + + +ippStage +parseStage(psString stageStr) +{ + if (!strcmp("raw", stageStr)) { + return IPP_STAGE_RAW; + } else if (!strcmp("chip", stageStr)) { + return IPP_STAGE_CHIP; + } else if (!strcmp("warp", stageStr)) { + return IPP_STAGE_WARP; + } else if (!strcmp("diff", stageStr)) { + return IPP_STAGE_DIFF; + } else { + psError(PS_ERR_UNKNOWN, true, "unknown stage string: %s\n", stageStr); + bail_out("", PS_EXIT_DATA_ERROR); + // notreached + return IPP_STAGE_NONE; + } +} + + +pmConfig *parseArguments(int argc, char **argv) { + + pmConfig *config = pmConfigRead(&argc, argv, NULL); + if (!config) { + bail_out("failed to parse arguments", PS_EXIT_CONFIG_ERROR); + } + + // bool status; + int argnum; + ippStage stage = IPP_STAGE_NONE; + if ((argnum = psArgumentGet(argc, argv, "-stage"))) { + psArgumentRemove(argnum, &argc, argv); + stage = parseStage(argv[argnum]); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "STAGE", 0, + "pipeline stage for streak removal", stage); + psArgumentRemove(argnum, &argc, argv); + } else { + psError(PS_ERR_UNKNOWN, true, "-stage is required\n"); + return NULL; + } + + if ((argnum = psArgumentGet(argc, argv, "-class_id"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "CLASS_ID", 0, + "class_id", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else if ((stage == IPP_STAGE_RAW) || (stage == IPP_STAGE_CHIP)) { + psError(PS_ERR_UNKNOWN, true, "-class_id is required\n"); + return NULL; + } + + if ((argnum = psArgumentGet(argc, argv, "-streaks"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STREAKS", 0, + "STREAKS", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else { + psError(PS_ERR_UNKNOWN, true, "-streaks is required\n"); + return NULL; + } + + if ((argnum = psArgumentGet(argc, argv, "-image"))) { + // for raw and chip level images we use psFits for I/O and ... + if (CHIP_LEVEL_INPUT(stage)) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT", 0, + "name of input image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else { + // ... for warped images we use pmFPAfiles + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "INPUT", "-image", NULL)) { ; + psError(PS_ERR_UNKNOWN, false, "failed to process -image"); + return NULL; + } + } + } else { + psError(PS_ERR_UNKNOWN, true, "-image is required\n"); + return NULL; + } + + + if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", NULL)) { ; + if ((stage == IPP_STAGE_RAW) || (stage == IPP_STAGE_CHIP)) { + psError(PS_ERR_UNKNOWN, true, "-astrom is required for raw and chip stages\n"); + return NULL; + } + } + if ((argnum = psArgumentGet(argc, argv, "-mask"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "MASK", 0, + "name of input mask image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } + if ((argnum = psArgumentGet(argc, argv, "-weight"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "WEIGHT", 0, + "name of input weight image", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } + if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "root name of temporary files", + argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else { + psError(PS_ERR_UNKNOWN, true, "-outroot is required\n"); + return NULL; + } + if ((argnum = psArgumentGet(argc, argv, "-recovery"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "Name of recovery root", + argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } else if (stage == IPP_STAGE_RAW) { + psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw\n"); + return NULL; + } + + return config; +} + +static void +copyPHU(streakFiles *sfiles) +{ + psAssert(sfiles->stage == IPP_STAGE_RAW, "copyPHU should only be used for raw stage"); + + psMetadata *imageHeader = psFitsReadHeader(NULL, sfiles->inImage->fits); + if (!imageHeader) { + psError(PS_ERR_IO, false, "failed to read primary header from %s\n", + sfiles->inImage->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + + // TODO: add keyword indicating that streaks have been removed + if (!psFitsWriteBlank(sfiles->outImage->fits, imageHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->outImage->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + // TODO: add keyword indicating that this is the recovery image + if (!psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->recImage->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + psFree(imageHeader); + + sFile *inMask = sfiles->inMask; + if (inMask) { + psMetadata *maskHeader = psFitsReadHeader(NULL, inMask->fits); + if (!maskHeader) { + psError(PS_ERR_IO, false, "failed to read primary header from %s\n", + sfiles->inMask->resolved_name); + bail_out("", 1); + } + // TODO: add keyword indicating that streaks have been removed + if (!psFitsWriteBlank(sfiles->outMask->fits, maskHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->outMask->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + // TODO: add keyword indicating that this is the recovery image + if (!psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->recMask->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + psFree(maskHeader); + } + sFile *inWeight = sfiles->inWeight; + if (inWeight) { + psMetadata *weightHeader = psFitsReadHeader(NULL, inWeight->fits); + if (!weightHeader) { + psError(PS_ERR_IO, false, "failed to read primary header from %s\n", + sfiles->inWeight->resolved_name); + bail_out("", 1); + } + // TODO: add keyword indicating that streaks have been removed + if (!psFitsWriteBlank(sfiles->outWeight->fits, weightHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->outWeight->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + // TODO: add keyword indicating that this is a recovery image + if (!psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) { + psError(PS_ERR_IO, false, "failed to write primary header to %s", + sfiles->recWeight->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } + psFree(weightHeader); + } +} + +// Determine whether the current header for this file is an image or not +static bool +isImage(sFile *in) +{ + bool status; + psString exttype = psMetadataLookupStr(&status, in->header, "EXTTYPE"); + if (exttype && !strcmp(exttype, "IMAGE")) { + return true; + } + // examine current header and determine if it is an image + psString xtension = psMetadataLookupStr(&status, in->header, "XTENSION"); + if (xtension) { + if (!strcmp(xtension, "IMAGE")) { + return true; + } else if (!strcmp(xtension, "BINTABLE")) { + psString ztension = psMetadataLookupStr(&status, in->header, "ZTENSION"); + if (ztension && !strcmp(ztension, "IMAGE")) { + return true; + } + } + } else if (in->nHDU == 1) { + // no extensions in the file, can just return true? For now require + // at least one dimension + int naxis = psMetadataLookupS32(&status, in->header, "NAXIS"); + return naxis > 0; + } + return false; +} + +static void +readImageFrom_pmFile(streakFiles *sf) +{ + // XXX: Currently this function assumes that it is only used for a single + // chip single cell FPA (i.e. a skyfile) + pmFPAview *view = sf->view; + assert(view->chip == 0); + view->cell = 0; + sf->cell = pmFPAviewThisCell(view, sf->inImage->pmfile->fpa); + + if (!sf->cell) { + bail_out("no cells found in chip", PS_EXIT_PROG_ERROR); + } + if (sf->cell->readouts->n != 1) { + psError(PS_ERR_PROGRAMMING, true, "unexpected number of readouts found: %ld", sf->cell->readouts->n); + bail_out("", PS_EXIT_PROG_ERROR); + } + view->readout = 0; + pmReadout *readout = pmFPAviewThisReadout(view, sf->inImage->pmfile->fpa); + if (!readout) { + bail_out("readout 0 not found", PS_EXIT_PROG_ERROR); + } + + sf->inImage->image = (psImage*) psMemIncrRefCounter(readout->image); +} + +static void +readImage(sFile *in, int extnum) +{ + psRegion region = {0, 0, 0, 0}; + + if (in->header) psFree(in->header); + if (in->image) { + psFree(in->image); + in->image = NULL; + } + if (in->imagecube) { + psFree(in->imagecube); + in->imagecube = NULL; + } + + in->header = psFitsReadHeader(NULL, in->fits); + if (!in->header) { + psError(PS_ERR_IO, false, "failed to read header from %s extnum: %d", + in->resolved_name, extnum); + bail_out("", PS_EXIT_DATA_ERROR); + } + + + if (!isImage(in)) { + return; + } + + in->numZPlanes = psMetadataLookupS32(NULL, in->header, "NAXIS3"); + + if (in->numZPlanes == 0) { + in->image = psFitsReadImage(in->fits, region, 0); + if (!in->image) { + psError(PS_ERR_IO, false, "failed to read image from %s extnum: %d", + in->resolved_name, extnum); + bail_out("", PS_EXIT_DATA_ERROR); + } + in->numCols = in->image->numCols; + in->numRows = in->image->numRows; + } else { + in->imagecube = psFitsReadImageCube(in->fits, region); + if (!in->imagecube) { + psError(PS_ERR_IO, false, "failed to read image cube from %s extnum: %d", + in->resolved_name, extnum); + bail_out("", PS_EXIT_DATA_ERROR); + } + psImage *image = (psImage *) (in->imagecube->data[0]); + in->numCols = image->numCols; + in->numRows = image->numRows; + } + +} + +static void +setOptions(sFile *sfile, psString extname, int bitpix, float bscale, float bzero) +{ + if (sfile->fits->options) { + psFree(sfile->fits->options); + } + sfile->fits->options = psFitsOptionsAlloc(); + sfile->fits->options->extword = psStringCopy(extname); + sfile->fits->options->scaling = PS_FITS_SCALE_MANUAL; + sfile->fits->options->bitpix = bitpix; + sfile->fits->options->bscale = bscale; + sfile->fits->options->bzero = bzero; +} + +static void +copyFitsOptions(sFile *out, sFile *rec, sFile *in) +{ + // Get current BITPIX, BSCALE, BZERO, EXTNAME + // Probably not necessary to look the numerical values up in this + // way, but guards against changes to psLib and cfitsio FITS + // handling. + psMetadataItem *bitpixItem = psMetadataLookup(in->header, "BITPIX"); + psAssert(bitpixItem, "Every FITS image should have BITPIX"); + int bitpix = psMetadataItemParseS32(bitpixItem); + psMetadataItem *bscaleItem = psMetadataLookup(in->header, "BSCALE"); + + float bscale; + if (!bscaleItem) { + psWarning("BSCALE isn't set; defaulting to unity"); + bscale = 1.0; + } else { + bscale = psMetadataItemParseF32(bscaleItem); + } + psMetadataItem *bzeroItem = psMetadataLookup(in->header, "BZERO"); + float bzero; + if (!bzeroItem) { + psWarning("BZERO isn't set; defaulting to zero"); + bzero = 0.0; + } else { + bzero = psMetadataItemParseF32(bzeroItem); + } + bool mdok; + + psString extname = psMetadataLookupStr(&mdok, in->header, "EXTNAME"); + setOptions(out, extname, bitpix, bscale, bzero); + setOptions(rec, extname, bitpix, bscale, bzero); +} + + +static void +copyTable(sFile *out, sFile *in, int extnum) +{ + bool mdok; + psString xtension = psMetadataLookupStr(&mdok, in->header, "XTENSION"); + psString extname = psMetadataLookupStr(&mdok, in->header, "EXTNAME"); + + if (!xtension) { + psWarning("extnum %d has no image and xtension is not defined in %s", + extnum, in->resolved_name); + // XXX abort? + return; + } + if (!extname) { + psWarning("extnum %d has no image and extname not defined in %s", + extnum, in->resolved_name); + // XXX abort? + return; + } + psArray *table = psFitsReadTable(in->fits); + if (!table) { + psError(PS_ERR_UNKNOWN, false, "failed to read table in extension %d from in->resolved name", extnum); + bail_out("", PS_EXIT_DATA_ERROR); + } + + if (!psFitsWriteTable(out->fits, out->header, table, extname)) { + psError(PS_ERR_UNKNOWN, false, "failed to copy table in extension %d", extnum); + bail_out("", PS_EXIT_DATA_ERROR); + } +} + +static void +createRecoveryImage(sFile *rec, sFile *in, int extnum) +{ + // PFS: may not be necessary but replaced in->numRows with in->image->numRows + rec->image = psImageRecycle(rec->image, in->image->numCols, in->image->numRows, in->image->type.type); + if (!rec->image) { + psError(PS_ERR_UNKNOWN, false, "failed to allocate recoveyr image for extnsion %d", extnum); + bail_out("", PS_EXIT_UNKNOWN_ERROR); + } + psImageInit(rec->image, 0); +} + +// set the image for the output files to the contents of the corresponding input file. +static bool +readAndCopyToOutput(streakFiles *sf) +{ + if (sf->inImage->pmfile) { + // image data from pmFPAfile (diff or warp file) + readImageFrom_pmFile(sf); + sf->astrom = streakSetAstrometry(sf->astrom, sf->inAstrom->fpa, sf->chip, true, sf->cell->concepts, + sf->inImage->numCols, sf->inImage->numRows); + } else { + // image data directly from psFits + readImage(sf->inImage, sf->extnum); + if (SFILE_IS_IMAGE(sf->inImage)) { + sf->astrom = streakSetAstrometry(sf->astrom, sf->inAstrom->fpa, sf->chip, false, sf->inImage->header, + sf->inImage->numCols, sf->inImage->numRows); + } + } + sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); + sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); + + if (!SFILE_IS_IMAGE(sf->inImage)) { + // fprintf(stderr, "skipping extnum: %d not an image\n", sf->extnum); + // XXX: This extension is probably a video table. + copyTable(sf->outImage, sf->inImage, sf->extnum); + return false; + } + + if (sf->inImage->image) { + sf->outImage->image = (psImage*) psMemIncrRefCounter(sf->inImage->image); + createRecoveryImage(sf->recImage, sf->inImage, sf->extnum); + } else { + // we don't do a reference if we have an imagecube + } + + // set up the compression parameters + copyFitsOptions(sf->outImage, sf->recImage, sf->inImage); + + // XXX: TODO: can we derive these values from the input header? + // psFitsCompressionGet(sf->inImage->image) gives compression none + // perhaps we should just use the definition of COMP_IMG in the configuration + psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); + psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); + + if (sf->inMask) { + readImage(sf->inMask, sf->extnum); + sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); + sf->outMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image); + sf->recMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image); + + createRecoveryImage(sf->recMask, sf->inMask, sf->extnum); + copyFitsOptions(sf->outMask, sf->recMask, sf->inMask); + + // XXX: see note above + psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); + psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); + } + + if (sf->inWeight) { + readImage(sf->inWeight, sf->extnum); + sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); + sf->outWeight->image = (psImage*) psMemIncrRefCounter(sf->inWeight->image); + + sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); + createRecoveryImage(sf->recWeight, sf->inWeight, sf->extnum); + + copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight); + // XXX: see note above + psFitsSetCompression(sf->outWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); + psFitsSetCompression(sf->recWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0); + } + + // and for raw images, create sub images that represent the actual image + // area (no overscan) + + + return true; +} + +static void +writeImage(sFile *sfile, psString extname, int extnum) +{ + if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) { + psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", + sfile->resolved_name, extnum); + bail_out("", PS_EXIT_DATA_ERROR); + } + psFree(sfile->image); + sfile->image = NULL; + psFree(sfile->header); + sfile->header = NULL; +} +static void +writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum) +{ + if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) { + psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", + sfile->resolved_name, extnum); + bail_out("", PS_EXIT_DATA_ERROR); + } + psFree(sfile->header); + sfile->header = NULL; +} + +static void +writeImages(streakFiles *sf) +{ + psString extname = NULL; + if (sf->nHDU > 1) { + extname = psMetadataLookupStr(NULL, sf->inImage->header, "EXTNAME"); + } + if (sf->inImage->numZPlanes == 0) { + writeImage(sf->outImage, extname, sf->extnum); + writeImage(sf->recImage, extname, sf->extnum); + } else { + // XXX: TODO: if a streak touched this cell then + // write the image cube to the recovery image if no streaks write it + // to the output image + writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum); + } + if (sf->outMask) { + writeImage(sf->outMask, extname, sf->extnum); + writeImage(sf->recMask, extname, sf->extnum); + } + if (sf->outWeight) { + writeImage(sf->outWeight, extname, sf->extnum); + writeImage(sf->recWeight, extname, sf->extnum); + } +} + +static void +closeImage(sFile *sfile) +{ + if (!sfile->fits) { + return; + } + if (!psFitsClose(sfile->fits)) { + psError(PS_ERR_IO, false, "failed to close image to %s", + sfile->resolved_name); + bail_out("", PS_EXIT_DATA_ERROR); + } +} + +static void +closeImages(streakFiles *sf) +{ + closeImage(sf->inImage); + closeImage(sf->outImage); + closeImage(sf->recImage); + if (sf->inMask) { + closeImage(sf->inMask); + closeImage(sf->outMask); + closeImage(sf->recMask); + } + if (sf->inWeight) { + closeImage(sf->inWeight); + closeImage(sf->outWeight); + closeImage(sf->recWeight); + } +} + +void +dumpAstro(streakFiles *sf) +{ + strkPt pt; + if (!sf->inImage->image) { + return; + } + + cellToSky(&pt, sf->astrom, 0, 0); + + psString extname = psMetadataLookupStr(NULL, sf->inImage->header, "EXTNAME"); + + printf("%2d %s %f %f\n", sf->extnum, extname, pt.x, pt.y); +} + + +int +main(int argc, char *argv[]) +{ + long i; + double imageValue, weightValue, maskValue; + StreakPixels *pixels; + PixelPos *pixelPos; + pmConfig *config = parseArguments(argc, argv); + if (!config) { + psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n"); + return PS_EXIT_CONFIG_ERROR; + } + + psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS"); + Streaks *streaks = readStreaksFile(streaksFileName); + + streakFiles *sfiles = openFiles(config); + + if (sfiles->stage == IPP_STAGE_RAW) { + // copy PHU to output files + copyPHU(sfiles); + // advance to the first image extension + if (!streakFilesNextExtension(sfiles)) { + psErrorStackPrint(stderr, "failed to advance to first extension of input images"); + exit (PS_EXIT_PROG_ERROR); + } + } + + // Iterate through components of image + do { + + // read the images and copy the pixels image from the inputs to the outputs + // This also sets up the astrometry and initializes the recovery images + if (!readAndCopyToOutput(sfiles)) { + // this extension is not an image skip (video descriptor? + // XXX is there anything else that could be in an input file that we need to handle? + continue; + } + + // Identify pixels to mask because of streaks + + pixels = streak_on_component (streaks, sfiles->astrom, + sfiles->inImage->numCols, sfiles->inImage->numRows); + +#ifdef notyet + // XXX: How is this formatted? + // PFS: Aren't these chips that are never warped? They should be cleared + // entirely or copy them to archive location. I don't know how + // to find these chips. The fact that a warp is not created for + // a chip must be recorded somewhere. I don't like this approach + // of including everything (ALL) and then removing them if they + // are warped. I don't know if it is practical. + // From ICD: + // In the raw and detrended images, the pixels which were not + // included in any of the streak-processed warps must also be masked. + // Note that the warp and difference skycells are only generated if more + // than a small fraction of the pixels are lit by the input image. + + // Identify pixels to mask because not in a warp + warp_pixels = pixel_list_initialise(ALL); // List of pixels to be excised because not in a warp + foreach warp (warps) { + // Calculate warp region on image + image_warp = warp_on_component(warp, astrom, numCols, numRows); + foreach pixel specified by image_warp { + remove_pixel(warp_pixels, pixel); + } + } +#endif + + for (i = 0; i < psArrayLength (pixels); ++i) { + pixelPos = psArrayGet (pixels, i); + imageValue = psImageGet (sfiles->inImage->image, pixelPos->x, pixelPos->y); + weightValue = psImageGet (sfiles->inWeight->image, pixelPos->x, pixelPos->y); + maskValue = psImageGet (sfiles->inMask->image, pixelPos->x, pixelPos->y); + + psImageSet (sfiles->recImage->image, pixelPos->x, pixelPos->y, imageValue); + psImageSet (sfiles->recWeight->image, pixelPos->x, pixelPos->y, weightValue); + psImageSet (sfiles->recMask->image, pixelPos->x, pixelPos->y, maskValue); + + psImageSet (sfiles->outImage->image, pixelPos->x, pixelPos->y, NAN); + psImageSet (sfiles->outWeight->image, pixelPos->x, pixelPos->y, NAN); + // TODO: + // Need to get mask weight for PM_MASK_DETECTOR. Is this stored + // in a config somewhere? Not sure how to properly set the maskValue. + psImageSet (sfiles->outMask->image, pixelPos->x, pixelPos->y, maskValue); + } + psArrayElementsFree (pixels); + + // dumpAstro(sfiles); + // write out the destreaked temporary images + writeImages(sfiles); + } while (streakFilesNextExtension(sfiles)); + + // close the destreaked temporary images + closeImages(sfiles); + + // TODO: + // Now copy the destreaked images files that we made above + // over each nebulous instances of input files + + // NOTE: from here on we can't just quit if something goes wrong. + // especially if we're working on a rawImage + + + return 0; +}