Index: trunk/magic/remove/src/.cvsignore =================================================================== --- trunk/magic/remove/src/.cvsignore (revision 21155) +++ trunk/magic/remove/src/.cvsignore (revision 21156) @@ -2,2 +2,3 @@ streaksreplace streakscompare +streaksrelease Index: trunk/magic/remove/src/Line.c =================================================================== --- trunk/magic/remove/src/Line.c (revision 21155) +++ trunk/magic/remove/src/Line.c (revision 21156) @@ -321,9 +321,11 @@ if (DistanceSquared (line, x, y) <= halfWidth2) { - pixel = psAlloc (sizeof(PixelPos)); - pixel->x = (int) x; - pixel->y = (int) y; - psArrayAdd (pixels, 1024, pixel); - psFree (pixel); + if (x >=0 && y >= 0) { + pixel = psAlloc (sizeof(PixelPos)); + pixel->x = (unsigned int) x; + pixel->y = (unsigned int) y; + psArrayAdd (pixels, 1024, pixel); + psFree (pixel); + } } } @@ -368,9 +370,11 @@ if (DistanceSquared (line, x, y) <= halfWidth2) { - pixel = psAlloc (sizeof(PixelPos)); - pixel->x = (int) x; - pixel->y = (int) y; - psArrayAdd (pixels, 1024, pixel); - psFree (pixel); + if (x >=0 && y >= 0) { + pixel = psAlloc (sizeof(PixelPos)); + pixel->x = (unsigned int) x; + pixel->y = (unsigned int) y; + psArrayAdd (pixels, 1024, pixel); + psFree (pixel); + } } } Index: trunk/magic/remove/src/Makefile.simple =================================================================== --- trunk/magic/remove/src/Makefile.simple (revision 21155) +++ trunk/magic/remove/src/Makefile.simple (revision 21156) @@ -6,10 +6,10 @@ COMMON_OBJECTS = \ streaksio.o \ - streaksutil.o + streaksutil.o \ + streaksastrom.o REMOVE_OBJECTS= \ ${COMMON_OBJECTS} \ streaksremove.o \ - streaksastrom.o \ streaksextern.o \ warpedpixels.o \ @@ -30,6 +30,6 @@ STREAKSFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1 OPTFLAGS= -g -#CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS} -CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} +CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS} +#CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} LDFLAGS=`psmodules-config --libs` Index: trunk/magic/remove/src/streaksastrom.c =================================================================== --- trunk/magic/remove/src/streaksastrom.c (revision 21155) +++ trunk/magic/remove/src/streaksastrom.c (revision 21156) @@ -17,5 +17,9 @@ if (fromCell) { - // The metadata is the cell concepts + // XXX: I think that this is not needed + // We are only asked to get the concepts from the cell when we're reading from a pmfile (skycell) + // In that case we (x0,y0) = (0,0) and the parirties are positive + // If the x0, y0 are not zero they are probably in the FITS unity based array format + // We should probably assert that these are zero cell_x0 = psMetadataLookupS32(&mdok, md, "CELL.X0"); if (!mdok) { @@ -28,5 +32,4 @@ return NULL; } - xParityCell = psMetadataLookupS32(&mdok, md, "CELL.XPARITY"); if (!mdok || (xParityCell != 1 && xParityCell != -1)) { @@ -65,7 +68,7 @@ return false; } - } else { // no metadata regular cell + // used by warpedpixels cell_x0 = 0; cell_y0 = 0; @@ -89,4 +92,50 @@ } +void +chipToCell(double *xCell, double *yCell, strkAstrom *astrom, double xChip, double yChip) +{ + // convert from chip to cell + if (astrom->xParity > 0) { + *xCell = xChip - astrom->cell_x0; + } else { + *xCell = astrom->cell_x0 - 1 - xChip; + } + if (astrom->yParity > 0) { + *yCell = yChip - astrom->cell_y0; + } else { + *yCell = astrom->cell_y0 - 1 - yChip; + } +} + +void +cellToChip(double *xChip, double *yChip, strkAstrom *astrom, double xCell, double yCell) +{ + if (astrom->xParity > 0) { + *xChip = xCell + astrom->cell_x0; + } else { + *xChip = astrom->cell_x0 - 1 - xCell; + } + if (astrom->yParity > 0) { + *yChip = yCell + astrom->cell_y0; + } else { + *yChip = astrom->cell_y0 - 1 - yCell; + } +} + +void +cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell) +{ + if (astrom->xParity > 0) { + *xChip = xCell + astrom->cell_x0; + } else { + *xChip = astrom->cell_x0 - 1 - xCell; + } + if (astrom->yParity > 0) { + *yChip = yCell + astrom->cell_y0; + } else { + *yChip = astrom->cell_y0 - 1 - yCell; + } +} + bool skyToCell(strkPt *outPt, strkAstrom *astrom, double ra, double dec) @@ -107,7 +156,5 @@ 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; + chipToCell(&outPt->x, &outPt->y, astrom, pt->chip->x, pt->chip->y); // printf("cell: %f %f chip: %f %f\n", outPt->x, outPt->y, pt->chip->x, pt->chip->y); @@ -115,5 +162,5 @@ return true; } - + bool cellToSky(strkPt *outPt, strkAstrom *astrom, double x, double y) @@ -123,7 +170,5 @@ 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; + cellToChip(&pt->chip->x, &pt->chip->y, astrom, x, y); pt->chip->xErr = 0; // Is setting these errors to zero the right thing to do? @@ -142,11 +187,4 @@ } -void -cellToChip(PixelPos *chip, strkAstrom *astrom, PixelPos *cell) -{ - // TODO: do I need to worry about going from int to float and back again here? - chip->x = (cell->x * astrom->xParity) + astrom->cell_x0; - chip->y = (cell->y * astrom->yParity) + astrom->cell_y0; -} static bool @@ -166,5 +204,4 @@ PS_ASSERT_PTR_NON_NULL(hdu->header, 1) if (sf->bilevelAstrometry) { - // Do we get here for GPC1 ? I don't necessarily have an fpa for the input image if (!pmAstromReadBilevelMosaic(sf->inAstrom->fpa, phu->header)) { psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA."); @@ -232,2 +269,10 @@ } } + +void +linearizeTransforms(strkAstrom *astrom) +{ + if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip)) { + streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR); + } +} Index: trunk/magic/remove/src/streakscompare.c =================================================================== --- trunk/magic/remove/src/streakscompare.c (revision 21155) +++ trunk/magic/remove/src/streakscompare.c (revision 21156) @@ -39,6 +39,6 @@ streaksExit("failed to advance to next extesion\n", PS_EXIT_DATA_ERROR); } - readImage(file1, 0, stage); - readImage(file2, 0, stage); + readImage(file1, 0, stage, false); + readImage(file2, 0, stage, false); psImage *image1 = file1->image; Index: trunk/magic/remove/src/streaksio.c =================================================================== --- trunk/magic/remove/src/streaksio.c (revision 21155) +++ trunk/magic/remove/src/streaksio.c (revision 21156) @@ -41,5 +41,5 @@ sf->inMask = sFileOpen(config, stage, "INPUT.MASK", NULL, false); - if (sf->inMask) { + if (sf->inMask && (sf->stage != IPP_STAGE_RAW)) { inputBasename = basename(sf->inMask->name); sf->outMask = sFileOpen(config, stage, "OUTPUT", inputBasename, true); @@ -295,5 +295,5 @@ streakFilesNextExtension(streakFiles *sf) { - // unless stage is raw or chip when we get here we're done + // unless stage is raw when we get here we're done if (sf->stage != IPP_STAGE_RAW) { if (sf->view) { @@ -315,12 +315,8 @@ 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); - } + // No mask and weight images for raw stage return true; } else { + // we're all done return false; } @@ -354,5 +350,5 @@ sFile *inMask = sfiles->inMask; - if (inMask) { + if (inMask && sfiles->outMask) { psMetadata *maskHeader = psFitsReadHeader(NULL, inMask->fits); if (!maskHeader) { @@ -469,8 +465,8 @@ } -void -setDataExtent(ippStage stage, sFile *in) -{ - if (stage == IPP_STAGE_RAW) { +static void +setDataExtent(ippStage stage, sFile *in, bool rawImage) +{ + if (rawImage) { psString datasec = psMetadataLookupStr(NULL, in->header, "DATASEC"); if (!datasec) { @@ -480,4 +476,19 @@ int xmin, xmax, ymin, ymax; sscanf(datasec, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax); + + // I've inadvertantly introduced an implicit assumption that cell coordinates + // are the same as image coordinates. This is true for GPC1 + // Test for this condition not being true + if (xmin != 1) { + psError(PS_ERR_IO, true, "can't handle datasec with minimum x != 1: %s", datasec); + streaksExit("", PS_EXIT_PROG_ERROR); + } + if (ymin != 1) { + psError(PS_ERR_IO, true, "can't handle datasec with minimum y != 1: %s", datasec); + streaksExit("", PS_EXIT_PROG_ERROR); + } + + // These values are different than the size of the actual image which includes + // the bias area. in->numCols = xmax - xmin + 1; in->numRows = ymax - ymin + 1; @@ -489,5 +500,5 @@ void -readImage(sFile *in, int extnum, ippStage stage) +readImage(sFile *in, int extnum, ippStage stage, bool isMask) { psRegion region = {0, 0, 0, 0}; @@ -526,4 +537,11 @@ streaksExit("", PS_EXIT_DATA_ERROR); } + if (in->image->type.type == PS_TYPE_U16) { + in->exciseValue = 65535; + psMetadataAddU16(in->header, PS_LIST_TAIL, "BLANK", 0, "", 65535); + psMetadataAddU16(in->header, PS_LIST_TAIL, "ZBLANK", 0, "", 65535); + } else { + in->exciseValue = NAN; + } } else { if (stage != IPP_STAGE_RAW) { @@ -539,5 +557,5 @@ psImage *image = (psImage *) (in->imagecube->data[0]); } - setDataExtent(stage, in); + setDataExtent(stage, in, (stage == IPP_STAGE_RAW) && !isMask); } @@ -676,5 +694,5 @@ writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum) { - if (!imagecube) { + if (!sfile || !imagecube) { return; } @@ -687,58 +705,4 @@ sfile->header = NULL; } - -#ifdef notdef -void -writeImages(streakFiles *sf, bool exciseImageCube) -{ - psString extname = NULL; - if (sf->nHDU > 1) { - bool mdok; - extname = psMetadataLookupStr(&mdok, sf->inImage->header, "EXTNAME"); - } - if (sf->inImage->numZPlanes == 0) { - // note exciseing complete images is handled in readAndCopyToOutput - writeImage(sf->outImage, extname, sf->extnum); - writeImage(sf->recImage, extname, sf->extnum); - } else { - // we have an image cube - double initValue; - if (exciseImageCube) { - // copy the entire input image to the recovery image - writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum); - initValue = NAN; - } else { - // otherwise write it to the output - writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum); - initValue = 0; - } - - // borrow one of the images from the imagecube and set it to init value - psImage *image = psArrayGet (sf->inImage->imagecube, 0); - psMemIncrRefCounter(image); - psImageInit(image, initValue); - if (exciseImageCube) { - sf->outImage->image = image; - writeImage(sf->outImage, extname, sf->extnum); - } else { - // write zero valued image to reccovery - if (sf->recImage) { - 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) { - 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); - } -} -#endif void @@ -952,2 +916,71 @@ } +bool +isExciseValue(double value, double exciseValue) +{ + if (isnan(exciseValue)) { + return isnan(value); + } else { + return value == exciseValue; + } +} + +void +setMaskedToNAN(streakFiles *sfiles, psU8 maskMask, bool printCounts) +{ + int maskedPixels = 0; + int nandPixels = 0; + int nandWeights = 0; + + psImage *image = sfiles->outImage->image; + psImage *mask = sfiles->inMask->image; + psImage *weight = NULL; + if (sfiles->outWeight) { + weight = sfiles->outWeight->image; + } + double exciseValue = sfiles->inImage->exciseValue; + + if (printCounts) { + psTimerStart("NAN_MASKED"); + } + + for (int y=0; y < sfiles->inImage->numRows; y++) { + for (int x=0; x < sfiles->inImage->numCols; x++) { + // these gets are not necessary, we could just set the pixels to nan + // but I want to get the counts + double imageVal = psImageGet(image, x, y); + psU8 maskVal; + if (sfiles->stage == IPP_STAGE_RAW) { + int xChip, yChip; + cellToChipInt(&xChip, &yChip, sfiles->astrom, x, y); + maskVal = psImageGet(mask, xChip, yChip); + } else { + maskVal = psImageGet(mask, x, y); + } + if (maskVal & maskMask) { + ++maskedPixels; + if (!isExciseValue(imageVal, sfiles->inImage->exciseValue)) { + ++nandPixels; + psImageSet(image, x, y, exciseValue); + } + if (weight) { + double weightVal = weight ? psImageGet(weight, x, y) : 0; + if (!isnan(weightVal)) { + ++nandWeights; + psImageSet(weight, x, y, NAN); + } + } + } + } + } + if (printCounts) { + printf("time to NAN mask pixels: %f\n", psTimerClear("NAN_MASKED")); + int totalPixels = image->numRows * image->numCols; + printf("pixels: %10ld\n", totalPixels); + printf("masked pixels: %10ld %4.2f%%\n", maskedPixels, 100. * maskedPixels / totalPixels); + printf("nand pixels: %10ld %4.2f%%\n", nandPixels, 100. * nandPixels / totalPixels); + if (weight) { + printf("nand weights: %10ld %4.2f%%\n", nandWeights, 100. * nandWeights / totalPixels); + } + } +} Index: trunk/magic/remove/src/streaksio.h =================================================================== --- trunk/magic/remove/src/streaksio.h (revision 21155) +++ trunk/magic/remove/src/streaksio.h (revision 21156) @@ -12,5 +12,5 @@ void closeImage(sFile *sfile); -void readImage(sFile *sfile, int extnum, ippStage stage); +void readImage(sFile *sfile, int extnum, ippStage stage, bool isMask); void copyPHU(streakFiles *sfiles, bool remove); void copyTable(sFile *out, sFile *in, int extnum); @@ -25,3 +25,5 @@ bool streakFilesNextExtension(streakFiles *sf); +bool isExciseValue(double, double); +void setMaskedToNAN(streakFiles *sfiles, psU8 maskMask, bool printCounts); #endif // STREAKS_IO_H Index: trunk/magic/remove/src/streaksrelease.c =================================================================== --- trunk/magic/remove/src/streaksrelease.c (revision 21155) +++ trunk/magic/remove/src/streaksrelease.c (revision 21156) @@ -51,9 +51,4 @@ } - int totalPixels = 0; - int totalMaskedPixels = 0; - int totalNANdPixels = 0; - int totalNANdWeights = 0; - // Iterate through each component of the input (there is only one except for raw images) do { @@ -70,36 +65,5 @@ } - int maskedPixels = 0; - int nandPixels = 0; - int nandWeights = 0; - - psImage *image = sfiles->outImage->image; - psImage *mask = sfiles->outMask->image; - psImage *weight= sfiles->outWeight->image; - - for (int y=0; y < image->numRows; y++) { - for (int x=0; x < image->numCols; x++) { - // these gets are not necessary, we could just set the pixels to nan - // but I want to get the counts - double imageVal = psImageGet(image, x, y); - double weightVal = psImageGet(weight, x, y); - psU8 maskVal = psImageGet(mask, x, y); - if (maskVal & maskMask) { - ++maskedPixels; - if (!isnan(imageVal)) { - ++nandPixels; - psImageSet(image, x, y, NAN); - } - if (!isnan(weightVal)) { - ++nandWeights; - psImageSet(weight, x, y, NAN); - } - } - } - } - totalPixels += sfiles->inImage->numRows * sfiles->inImage->numCols; - totalMaskedPixels += maskedPixels; - totalNANdPixels += nandPixels; - totalNANdWeights += nandWeights; + setMaskedToNAN(sfiles, maskMask, true); // write out the destreaked temporary images and the recovery images @@ -108,9 +72,4 @@ printf("time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT")); } while (streakFilesNextExtension(sfiles)); - - printf("pixels: %10ld\n", totalPixels); - printf("masked pixels: %10ld %4.2f%%\n", totalMaskedPixels, 100. * totalMaskedPixels / totalPixels); - printf("nand pixels: %10ld %4.2f%%\n", totalNANdPixels, 100. * totalNANdPixels / totalPixels); - printf("nand weights: %10ld %4.2f%%\n", totalNANdWeights, 100. * totalNANdWeights / totalPixels); psTimerStart("CLOSE_IMAGES"); @@ -292,5 +251,5 @@ } else { // image data directly from psFits - readImage(sf->inImage, sf->extnum, sf->stage); + readImage(sf->inImage, sf->extnum, sf->stage, false); } sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); @@ -325,5 +284,5 @@ if (sf->inMask) { - readImage(sf->inMask, sf->extnum, sf->stage); + readImage(sf->inMask, sf->extnum, sf->stage, true); sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); if (sf->recMask) { @@ -358,5 +317,5 @@ if (sf->inWeight) { - readImage(sf->inWeight, sf->extnum, sf->stage); + readImage(sf->inWeight, sf->extnum, sf->stage, false); sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); if (sf->recWeight) { Index: trunk/magic/remove/src/streaksremove.c =================================================================== --- trunk/magic/remove/src/streaksremove.c (revision 21155) +++ trunk/magic/remove/src/streaksremove.c (revision 21156) @@ -9,4 +9,5 @@ static void writeImages(streakFiles *sf, bool exciseImageCube); static bool replicateOutputs(streakFiles *sfiles); +static void updateAstrometry(streakFiles *sfiles); int @@ -21,5 +22,5 @@ if (!config) { psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n"); - return PS_EXIT_CONFIG_ERROR; + streaksExit("", PS_EXIT_CONFIG_ERROR); } @@ -27,11 +28,19 @@ if (!status) { psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n"); - return PS_EXIT_CONFIG_ERROR; + streaksExit("", PS_EXIT_CONFIG_ERROR); } double maskStreak = (double) psMetadataLookupU8(&status, masks, "STREAK"); if (!status) { psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n"); - return PS_EXIT_CONFIG_ERROR; - } + streaksExit("", PS_EXIT_CONFIG_ERROR); + } + + // optionally setting pixels with any mask bits execpt POOR.WARP to NAN + psU8 poorWarp = (double) psMetadataLookupU8(&status, masks, "POOR.WARP"); + if (!status) { + psError(PM_ERR_CONFIG, false, "failed to lookup mask value for POOR.WARP in recipes\n"); + streaksExit("", PS_EXIT_CONFIG_ERROR); + } + psU8 maskMask = ~poorWarp; psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS"); @@ -40,9 +49,15 @@ if (!streaks) { psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName); - exit (PS_EXIT_PROG_ERROR); + streaksExit("", PS_EXIT_PROG_ERROR); } streakFiles *sfiles = openFiles(config, true); setupAstrometry(sfiles); + + bool nanForRelease = psMetadataLookupBool(&status, config->arguments, "NAN_FOR_RELEASE"); + if (nanForRelease && (sfiles->inMask == NULL)) { + psError(PS_ERR_UNKNOWN, true, "mask image must be supplied with -release."); + streaksExit("", PS_EXIT_PROG_ERROR); + } bool exciseAll = false; @@ -113,13 +128,13 @@ printf("time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS")); - if (checkNonWarpedPixels) { - psTimerStart("EXCISE_NON_WARPED"); - - exciseNonWarpedPixels(sfiles, maskStreak); - - printf("time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED")); - } if (sfiles->inImage->image) { + if (checkNonWarpedPixels) { + psTimerStart("EXCISE_NON_WARPED"); + + exciseNonWarpedPixels(sfiles, maskStreak); + + printf("time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED")); + } totalStreakPixels += psArrayLength(pixels); psTimerStart("REMOVE_STREAKS"); @@ -137,4 +152,9 @@ } printf("time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS")); + + if (nanForRelease) { + setMaskedToNAN(sfiles, maskMask, true); + } + } else { // this component contains an image cube, excise it completely @@ -142,4 +162,9 @@ } psArrayElementsFree (pixels); + } + + if (sfiles->stage == IPP_STAGE_CHIP) { + // if (CHIP_LEVEL_INPUT(sfiles->stage)) { + updateAstrometry(sfiles); } @@ -302,4 +327,10 @@ true); } + + if ((argnum = psArgumentGet(argc, argv, "-release"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "NAN_FOR_RELEASE", 0, "set masked pixels to NAN", + true); + } if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) { @@ -424,4 +455,23 @@ return config; +} + +static void +updateAstrometry(streakFiles *sf) +{ + if (sf->bilevelAstrometry) { + linearizeTransforms(sf->astrom); + // XXX: pmAstrometry.c should set the correct CTYPE values + if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) { + psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum); + streaksExit("", PS_EXIT_UNKNOWN_ERROR); + } + if (sf->outMask) { + pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); + } + if (sf->outWeight) { + pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001); + } + } } @@ -438,5 +488,5 @@ } else { // image data directly from psFits - readImage(sf->inImage, sf->extnum, sf->stage); + readImage(sf->inImage, sf->extnum, sf->stage, false); if (SFILE_IS_IMAGE(sf->inImage)) { sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, @@ -493,41 +543,43 @@ if (sf->inMask) { - readImage(sf->inMask, sf->extnum, sf->stage); - sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); - if (sf->recMask) { - sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); - } - if (updateAstrometry) { - pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); - } - setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll); - if (sf->outChMask) { - sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header); - sf->outChMask->image = (psImage *) psMemIncrRefCounter(sf->outMask->image); - if (sf->recChMask) { - sf->recChMask->header = (psMetadata *) psMemIncrRefCounter(sf->recMask->header); - sf->recChMask->image = (psImage *) psMemIncrRefCounter(sf->recMask->image); - } - } + readImage(sf->inMask, sf->extnum, sf->stage, true); + if (sf->outMask) { + sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); + if (sf->recMask) { + sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); + } + if (updateAstrometry) { + pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001); + } + setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll); + if (sf->outChMask) { + sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header); + sf->outChMask->image = (psImage *) psMemIncrRefCounter(sf->outMask->image); + if (sf->recChMask) { + sf->recChMask->header = (psMetadata *) psMemIncrRefCounter(sf->recMask->header); + sf->recChMask->image = (psImage *) psMemIncrRefCounter(sf->recMask->image); + } + } #ifdef STREAKS_COMPRESS_OUTPUT - // XXX: see note above - copyFitsOptions(sf->outMask, sf->recMask, sf->inMask); - psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - if (sf->recMask) { - psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - } - if (sf->outChMask) { - copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask); - psFitsSetCompression(sf->outChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - if (sf->recChMask) { - psFitsSetCompression(sf->recChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); - } - } + // XXX: see note above + copyFitsOptions(sf->outMask, sf->recMask, sf->inMask); + psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); + if (sf->recMask) { + psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); + } + if (sf->outChMask) { + copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask); + psFitsSetCompression(sf->outChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); + if (sf->recChMask) { + psFitsSetCompression(sf->recChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0); + } + } #endif + } } if (sf->inWeight) { - readImage(sf->inWeight, sf->extnum, sf->stage); + readImage(sf->inWeight, sf->extnum, sf->stage, false); sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); if (sf->recWeight) { @@ -661,4 +713,6 @@ excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue) { + double exciseValue = sfiles->inImage->exciseValue; + // we clip so that the streak calculation code doesn't have to if ((x < 0) || (x >= sfiles->inImage->numCols) || (y < 0) || (y >= sfiles->inImage->numRows)) { @@ -667,10 +721,10 @@ double imageValue = psImageGet (sfiles->inImage->image, x, y); - if (sfiles->recImage && !isnan(imageValue) ) { + if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) { psImageSet (sfiles->recImage->image, x, y, imageValue); } if (sfiles->transparentStreaks == 0) { - psImageSet (sfiles->outImage->image, x, y, NAN); + psImageSet (sfiles->outImage->image, x, y, exciseValue); } else { if (streak) { @@ -678,16 +732,17 @@ psImageSet (sfiles->outImage->image, x, y, imageValue + sfiles->transparentStreaks); } else { - psImageSet (sfiles->outImage->image, x, y, NAN); - } - } - - if (sfiles->inWeight) { + psImageSet (sfiles->outImage->image, x, y, exciseValue); + } + } + + if (sfiles->outWeight) { if (sfiles->recWeight) { double weightValue = psImageGet (sfiles->inWeight->image, x, y); psImageSet (sfiles->recWeight->image, x, y, weightValue); } + // Assume that weight images are always a floating point type psImageSet (sfiles->outWeight->image, x, y, NAN); } - if (sfiles->inMask) { + if (sfiles->outMask) { if (sfiles->recMask) { double maskValue = psImageGet (sfiles->inMask->image, x, y); @@ -705,9 +760,12 @@ int xParity = sfiles->astrom->xParity; int yParity = sfiles->astrom->yParity; - int numCols = sfiles->inImage->numCols; - int numRows = sfiles->inImage->numRows; + int numCols = sfiles->inImage->numCols; // for raw images this was calculated from the width of datasec + int numRows = sfiles->inImage->numRows; // for raw images this was calculated from the height of datasec // printf("%2d x0: %4d y0: %4d xpar: %d ypar: %d\n", sfiles->extnum, cell_x0, cell_y0, xParity, yParity); + // XXX: in the case of raw images we are making an implicit assumption here + // that the DATASEC begins at (x, y) = (0, 0) in the array for this image + // While this is true for GPC1 it isn't necessarily true for (int yCell = 0; yCell < numRows; yCell++) { int yChip; @@ -751,4 +809,5 @@ // we clip so that the streak calculation code doesn't have to + // clipping here insures that we don't touch the overscan regions if ((cellCoord->x < 0) || (cellCoord->x >= sfiles->inImage->numCols) || (cellCoord->y < 0) || (cellCoord->y >= sfiles->inImage->numRows)) { @@ -756,11 +815,6 @@ } - cellToChip(&chipCoord, sfiles->astrom, cellCoord); - -#ifdef notdef - // Shouldn't need this but for raw images I'm getting - // x = warpedPixels->numCols for cell xy70 when cell.x = 0 - // Perhaps cell_x0 should be IMNPIX-1 since IMNPIX is in the 1 based system of the fits - // headers + cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, cellCoord->x, cellCoord->y); + if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) { return false; @@ -769,5 +823,4 @@ return false; } -#endif return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false; Index: trunk/magic/remove/src/streaksremove.h =================================================================== --- trunk/magic/remove/src/streaksremove.h (revision 21155) +++ trunk/magic/remove/src/streaksremove.h (revision 21156) @@ -26,6 +26,11 @@ psArray *imagecube; pmFPAfile *pmfile; - int numCols; - int numRows; + int col0; + int numCols; // number of columns in data section of image (doesn't include bias section) + int xParity; + int row0; + int numRows; // number of rows in data section of image (does not include bias section) + int yParity; + double exciseValue; } sFile; @@ -75,5 +80,6 @@ int numCols, int numRows); // can't declare this in streaksastrom due to header file ordering -extern void cellToChip(PixelPos *chip, strkAstrom *astrom, PixelPos *cell); +extern void cellToChip(double *xChip, double *yChip, strkAstrom *astrom, double xCell, double yCell); +extern void cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell); extern bool computeWarpedPixels(streakFiles *sf); Index: trunk/magic/remove/src/streaksreplace.c =================================================================== --- trunk/magic/remove/src/streaksreplace.c (revision 21155) +++ trunk/magic/remove/src/streaksreplace.c (revision 21156) @@ -327,8 +327,8 @@ } else { // image data directly from psFits - readImage(sf->inImage, sf->extnum, sf->stage); + readImage(sf->inImage, sf->extnum, sf->stage, false); } // read the recovery pixels - readImage(sf->recImage, sf->extnum, sf->stage); + readImage(sf->recImage, sf->extnum, sf->stage, false); sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header); @@ -359,6 +359,6 @@ #endif if (sf->inMask) { - readImage(sf->inMask, sf->extnum, sf->stage); - readImage(sf->recMask, sf->extnum, sf->stage); + readImage(sf->inMask, sf->extnum, sf->stage, true); + readImage(sf->recMask, sf->extnum, sf->stage, true); sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header); @@ -372,6 +372,6 @@ if (sf->inWeight) { - readImage(sf->inWeight, sf->extnum, sf->stage); - readImage(sf->recWeight, sf->extnum, sf->stage); + readImage(sf->inWeight, sf->extnum, sf->stage, false); + readImage(sf->recWeight, sf->extnum, sf->stage, false); sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header); Index: trunk/magic/remove/src/warpedpixels.c =================================================================== --- trunk/magic/remove/src/warpedpixels.c (revision 21155) +++ trunk/magic/remove/src/warpedpixels.c (revision 21156) @@ -41,5 +41,7 @@ // write out the warped pixels psMetadata * header = psMetadataAlloc(); - pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001); + if (sf->inAstrom->fpa->toSky->type != PS_PROJ_DIS) { + pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001); + } psFits *fits = psFitsOpen("warpedpixels.fits", "w"); @@ -76,5 +78,4 @@ sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, NULL, sf->warpedPixels->numCols, sf->warpedPixels->numRows); - // convert corners of skycell to sky coordinates