Index: trunk/magic/remove/src/streaksremove.c =================================================================== --- trunk/magic/remove/src/streaksremove.c (revision 20520) +++ trunk/magic/remove/src/streaksremove.c (revision 20573) @@ -1,7 +1,7 @@ -#include "streaks.h" -#include "streaksextern.h" +#include "streaksremove.h" #include "libgen.h" #include "unistd.h" +// XXX: I don't think we need the lock and unlock functions extern bool sFileLock(sFile * sfile); extern bool sFileUnlock(sFile * sfile); @@ -10,5 +10,5 @@ static nebServer *ourNebServer = NULL; -// for clarity in this program we don't propagate errors up the stack +// to enhance clarity we don't propagate errors up the stack // we just bail out void streaksremoveExit(psString str, int exitCode) { @@ -70,5 +70,5 @@ } -psString +static psString resolveFilename(pmConfig *config, sFile *sfile, bool create) { @@ -166,10 +166,12 @@ sfile->resolved_name = resolveFilename(config, sfile, true); if (!sfile->resolved_name) { - psError(PS_ERR_IO, false, "Failed to create file %s", sfile->name); + psError(PS_ERR_IO, false, "Failed to resolve filename for %s", sfile->name); sFileFree(sfile); streaksremoveExit("", 1); } sfile->fits = psFitsOpen(sfile->resolved_name, "w"); - sfile->fits->options = psFitsOptionsAlloc(); + if (sfile->fits) { + sfile->fits->options = psFitsOptionsAlloc(); + } } else { sfile->name = psStringCopy(name); @@ -197,10 +199,11 @@ -bool -astrometry_read(streakFiles *sf) +static bool +readAstrometry(streakFiles *sf) { pmHDU *phu = pmFPAviewThisPHU(sf->view, sf->inAstrom->fpa); if (phu) { - char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + bool status; + char *ctype = psMetadataLookupStr(&status, phu->header, "CTYPE1"); if (ctype) { sf->bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); @@ -229,4 +232,5 @@ } +// load the fpa containing the astrometry, find the appropriate chip and process the data static void setupAstromFromFPA(streakFiles *sf) @@ -259,6 +263,5 @@ } - // read in the astrometry - astrometry_read(sf); + readAstrometry(sf); } @@ -312,7 +315,4 @@ } setupAstromFromFPA(sf); - if (CHIP_LEVEL_INPUT(sf->stage)) { - computeWarpedPixels(sf); - } psElemType tileType; // Type corresponding to "long" @@ -373,5 +373,5 @@ { // unless stage is raw or chip when we get here we're done - if (!CHIP_LEVEL_INPUT(sf->stage)) { + if (sf->stage != IPP_STAGE_RAW) { sf->view->readout = -1; pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER); @@ -442,5 +442,4 @@ } - // bool status; int argnum; ippStage stage = IPP_STAGE_NONE; @@ -667,4 +666,5 @@ // Determine whether the current header for this file is an image or not +// Find a cleaner way to do this static bool isImage(sFile *in) @@ -686,4 +686,6 @@ } } + } else if (psMetadataLookupBool(&status, in->header, "ZIMAGE")) { + return true; } else if (in->nHDU == 1) { // no extensions in the file, can just return true? For now require @@ -721,6 +723,25 @@ } +void +setDataExtent(ippStage stage, sFile *in) +{ + if (stage == IPP_STAGE_RAW) { + psString datasec = psMetadataLookupStr(NULL, in->header, "DATASEC"); + if (!datasec) { + psError(PS_ERR_IO, false, "failed to find DATASEC in header"); + streaksremoveExit("", PS_EXIT_DATA_ERROR); + } + int xmin, xmax, ymin, ymax; + sscanf(datasec, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax); + in->numCols = xmax - xmin + 1; + in->numRows = ymax - ymin + 1; + } else { + in->numCols = in->image->numCols; + in->numRows = in->image->numRows; + } +} + static void -readImage(sFile *in, int extnum) +readImage(sFile *in, int extnum, ippStage stage) { psRegion region = {0, 0, 0, 0}; @@ -756,7 +777,10 @@ streaksremoveExit("", PS_EXIT_DATA_ERROR); } - in->numCols = in->image->numCols; in->numRows = in->image->numRows; } else { + if (stage != IPP_STAGE_RAW) { + psError(PS_ERR_PROGRAMMING, true, "unexpected image cube found in non-raw image"); + streaksremoveExit("", PS_EXIT_PROG_ERROR); + } in->imagecube = psFitsReadImageCube(in->fits, region); if (!in->imagecube) { @@ -766,12 +790,10 @@ } psImage *image = (psImage *) (in->imagecube->data[0]); - in->numCols = image->numCols; - in->numRows = image->numRows; - } - + } + setDataExtent(stage, in); } static void -setFitsOptions(sFile *sfile, psString extname, int bitpix, float bscale, float bzero) +setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero) { if (sfile->fits->options) { @@ -812,9 +834,7 @@ bzero = psMetadataItemParseF32(bzeroItem); } - bool mdok; - - psString extname = psMetadataLookupStr(&mdok, in->header, "EXTNAME"); - setFitsOptions(out, extname, bitpix, bscale, bzero); - setFitsOptions(rec, extname, bitpix, bscale, bzero); + +// setFitsOptions(out, bitpix, bscale, bzero); +// setFitsOptions(rec, bitpix, bscale, bzero); } @@ -869,13 +889,17 @@ // 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); + sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, 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); + readImage(sf->inImage, sf->extnum, sf->stage); if (SFILE_IS_IMAGE(sf->inImage)) { - sf->astrom = streakSetAstrometry(sf->astrom, sf->inAstrom->fpa, sf->chip, false, + sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, sf->inImage->header, sf->inImage->numCols, sf->inImage->numRows); } + } + if (!sf->astrom) { + psError(PS_ERR_UNKNOWN, false, "failed to set up astrometry for extnsion %d", sf->extnum); + streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR); } sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); @@ -906,5 +930,5 @@ if (sf->inMask) { - readImage(sf->inMask, sf->extnum); + readImage(sf->inMask, sf->extnum, sf->stage); sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); sf->outMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image); @@ -920,5 +944,5 @@ if (sf->inWeight) { - readImage(sf->inWeight, sf->extnum); + readImage(sf->inWeight, sf->extnum, sf->stage); sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); sf->outWeight->image = (psImage*) psMemIncrRefCounter(sf->inWeight->image); @@ -943,4 +967,7 @@ writeImage(sFile *sfile, psString extname, int extnum) { + if (!sfile->image) { + return; + } if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) { psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", @@ -953,4 +980,5 @@ sfile->header = NULL; } + static void writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum) @@ -966,9 +994,10 @@ static void -writeImages(streakFiles *sf) +writeImages(streakFiles *sf, bool exciseAll) { psString extname = NULL; if (sf->nHDU > 1) { - extname = psMetadataLookupStr(NULL, sf->inImage->header, "EXTNAME"); + bool mdok; + extname = psMetadataLookupStr(&mdok, sf->inImage->header, "EXTNAME"); } if (sf->inImage->numZPlanes == 0) { @@ -976,8 +1005,26 @@ 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); + double initValue; + if (exciseAll) { + writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum); + initValue = NAN; + } else { + writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum); + initValue = 0; + } + + // borrow one of the images from the imagecube and set it to NAN + psImage *image = psArrayGet (sf->inImage->imagecube, 0); + psMemIncrRefCounter(image); + psImageInit(image, initValue); + if (exciseAll) { + sf->outImage->image = image; + writeImage(sf->outImage, extname, sf->extnum); + } else { + sf->recImage->image = image; + writeImage(sf->recImage, extname, sf->extnum); + } + // Assumption there are no mask and weight images for video cells + return; } if (sf->outMask) { @@ -1194,4 +1241,84 @@ } +static void +excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, double newMaskValue) +{ + double imageValue = psImageGet (sfiles->inImage->image, x, y); + psImageSet (sfiles->recImage->image, x, y, imageValue); + psImageSet (sfiles->outImage->image, x, 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. + + if (sfiles->inWeight) { + double weightValue = psImageGet (sfiles->inWeight->image, x, y); + psImageSet (sfiles->recWeight->image, x, y, weightValue); + psImageSet (sfiles->outWeight->image, x, y, NAN); + } + if (sfiles->inMask) { + double maskValue = psImageGet (sfiles->inMask->image, x, y); + psImageSet (sfiles->recMask->image, x, y, maskValue); + psImageSet (sfiles->outMask->image, x, y, newMaskValue); + } +} + +static void +exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue) +{ + int cell_x0 = sfiles->astrom->cell_x0; + int cell_y0 = sfiles->astrom->cell_y0; + int xParity = sfiles->astrom->xParity; + int yParity = sfiles->astrom->yParity; + int numCols = sfiles->inImage->numCols; + int numRows = sfiles->inImage->numRows; + +// printf("%2d x0: %4d y0: %4d xpar: %d ypar: %d\n", sfiles->extnum, cell_x0, cell_y0, xParity, yParity); + + for (int yCell = 0; yCell < numRows; yCell++) { + int yChip; + if (yParity == 1) { + yChip = cell_y0 + yCell; + } else { + yChip = cell_y0 - yCell; + } + + psU8 *pixels = sfiles->warpedPixels->data.U8[yChip]; + + if (xParity == 1) { + pixels += cell_x0; + for (int xCell = 0; xCell < numCols; xCell++, pixels++) { + if (! *pixels ) { + excisePixel(sfiles, xCell, yCell, newMaskValue); + } + } + } else { + // negative parity + // point to the right most pixel in chip coordinates (lowest in cell coords) + pixels += cell_x0 - numCols; + for (int xCell = numCols - 1; xCell >= 0 ; xCell--, pixels++) { + if (!*pixels) { + excisePixel(sfiles, xCell, yCell, newMaskValue); + } + } + } + } +} + +bool +warpedPixel(streakFiles *sfiles, PixelPos *cellCoord) +{ + PixelPos chipCoord; + + if (!CHIP_LEVEL_INPUT(sfiles->stage)) { + // if we're here on a skycell image by definition this pixel was warped + return true; + } + + cellToChip(&chipCoord, sfiles->astrom, cellCoord); + + return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false; +} + int main(int argc, char *argv[]) @@ -1220,38 +1347,5 @@ 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. + if (CHIP_LEVEL_INPUT(sfiles->stage)) { // From ICD: // In the raw and detrended images, the pixels which were not @@ -1260,39 +1354,67 @@ // 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); + computeWarpedPixels(sfiles); + } + + 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 { + bool exciseAll = false; + + // 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); + // XXX: + // + // PFS: Need to get mask weight for PM_MASK_DETECTOR. Is this stored + // in a config somewhere? Not sure how to properly set the maskValue. + // WES just pick some numbers +#define MASK_STREAK 42 +#define MASK_NOT_WARPED 51 + + exciseNonWarpedPixels(sfiles, MASK_NOT_WARPED); + + if (sfiles->inImage->image) { + for (i = 0; i < psArrayLength (pixels); ++i) { + pixelPos = psArrayGet (pixels, i); + + if (warpedPixel(sfiles, pixelPos)) { + + excisePixel(sfiles, pixelPos->x, pixelPos->y, MASK_STREAK); + + } else { + // This pixel was not included in any warp and has thus already excised above + } } - } -#endif - - // XXX: This loop following will SEGV for raw video Cells becuase the image is null - for (i = 0; i < psArrayLength (pixels); ++i) { - pixelPos = psArrayGet (pixels, i); - imageValue = psImageGet (sfiles->inImage->image, pixelPos->x, pixelPos->y); - psImageSet (sfiles->recImage->image, pixelPos->x, pixelPos->y, imageValue); - psImageSet (sfiles->outImage->image, pixelPos->x, pixelPos->y, NAN); - if (sfiles->inWeight) { - weightValue = psImageGet (sfiles->inWeight->image, pixelPos->x, pixelPos->y); - psImageSet (sfiles->recWeight->image, pixelPos->x, pixelPos->y, weightValue); - psImageSet (sfiles->outWeight->image, pixelPos->x, pixelPos->y, NAN); - } - if (sfiles->inMask) { - maskValue = psImageGet (sfiles->inMask->image, pixelPos->x, pixelPos->y); - psImageSet (sfiles->recMask->image, pixelPos->x, pixelPos->y, maskValue); - // 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); - } + } else { + // This video cell is touched by a streak excise + // For now we excise all video cells. + // if (psArrayLength(pixels)) + exciseAll = true; } psArrayElementsFree (pixels); // write out the destreaked temporary images and the recovery images - writeImages(sfiles); + writeImages(sfiles, exciseAll); + } while (streakFilesNextExtension(sfiles));