Index: trunk/magic/remove/src/Makefile.am =================================================================== --- trunk/magic/remove/src/Makefile.am (revision 26438) +++ trunk/magic/remove/src/Makefile.am (revision 26438) @@ -0,0 +1,60 @@ +if HAVE_SVNVERSION +STREAKS_VERSION=`$(SVNVERSION) ..` +else +STREAKS_VERSION="UNKNOWN" +endif + +if HAVE_SVN +STREAKS_BRANCH=`$(SVN) info .. | $(SED) -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'` +STREAKS_SOURCE=`$(SVN) info | $(SED) -n -e 's/Repository UUID: // p'` +else +STREAKS_BRANCH="UNKNOWN" +STREAKS_SOURCE="UNKNOWN" +endif + +# Force recompilation of streaksVersion.c, since it gets the version information +streaksVersion.c: streaksVersionDefinitions.h +streaksVersionDefinitions.h: streaksVersionDefinitions.h.in FORCE + -$(RM) streaksVersionDefinitions.h + $(SED) -e "s|@STREAKS_VERSION@|\"$(STREAKS_VERSION)\"|" -e "s|@STREAKS_BRANCH@|\"$(STREAKS_BRANCH)\"|" -e "s|@STREAKS_SOURCE@|\"$(STREAKS_SOURCE)\"|" streaksVersionDefinitions.h.in > streaksVersionDefinitions.h +FORCE: ; + +#bin_PROGRAMS = streaksremove streaksrelease streaksreplace streakscompare isdestreaked +bin_PROGRAMS = streaksremove + +STREAKSREMOVE_CFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1 + +streaksremove_CFLAGS = $(STREAKSREMOVE_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +streaksremove_LDFLAGS = $(STREAKSREMOVE_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) + +common_sources = \ + streaksio.c \ + streaksutil.c \ + streaksastrom.c + +streaksremove_SOURCES = \ + streaksremove.c \ + streaksextern.c \ + diffedpixels.c \ + Line.c \ + streaksio.c \ + streaksutil.c \ + streaksastrom.c \ + streaksVersion.c + +include_HEADERS = \ + streaksremove.h \ + streaksastrom.h \ + streaksextern.h \ + streaksio.h \ + streaksremove.h \ + streaksVersionDefinitions.h + +clean-local: + -rm -f TAGS + -rm -f tags + +# Tags for emacs and vi +tags: + etags `find . -name \*.[ch] -print` + ctags `find . -name \*.[ch] -print` Index: trunk/magic/remove/src/Makefile.simple =================================================================== --- trunk/magic/remove/src/Makefile.simple (revision 26408) +++ trunk/magic/remove/src/Makefile.simple (revision 26438) @@ -12,5 +12,5 @@ streaksremove.o \ streaksextern.o \ - warpedpixels.o \ + diffedpixels.o \ Line.o Index: trunk/magic/remove/src/diffedpixels.c =================================================================== --- trunk/magic/remove/src/diffedpixels.c (revision 26438) +++ trunk/magic/remove/src/diffedpixels.c (revision 26438) @@ -0,0 +1,332 @@ +#include "streaksremove.h" +#include "pmAstrometryWCS.h" + +// #define DEBUG_PRINT 1 + +static void addTouchedPixels(streakFiles *sf, psString fileName); +static void nameCorners(psPlane pt[4]); +static int xLeft(psPlane *pt, int y); +static int xRight(psPlane *pt, int y); + +// Examine the set of diff skycells and compute the pixels that were +// included in difference processing +bool +computeDiffedPixels(streakFiles *sf) +{ + bool status; + pmConfig *config = sf->config; + + psArray *skycells = psMetadataLookupPtr(&status, config->arguments, "SKYCELLS"); + if (skycells == NULL) { + return false; + } + + int n = psArrayLength(skycells); + + psRegion *bounds = pmChipPixels(sf->chip); + + int width = bounds->x1 - bounds->x0; + int height = bounds->y1 - bounds->y0; + + psFree(bounds); + + sf->diffedPixels = psImageAlloc(width, height, PS_TYPE_U8); + psImageInit(sf->diffedPixels, 0); + + for (int i=0; iinAstrom->fpa->toSky->type != PS_PROJ_DIS) { + pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001); + } + psFits *fits = psFitsOpen("diffedpixels.fits", "w"); + + psFitsWriteImage(fits, header, sf->diffedPixels, 0, NULL); + psFitsClose(fits); + psFree(header); + } + return true; +} + +static void +addTouchedPixels(streakFiles *sf, psString fileName) +{ + // read fits header + psFits *fits = psFitsOpen(fileName, "r"); + if (!fits) { + psError(PS_ERR_IO, false, "failed to open skycell file: %s", fileName); + streaksExit("", PS_EXIT_DATA_ERROR); + } + psMetadata *header = psFitsReadHeader(NULL, fits); + psFitsClose(fits); + if (!header) { + psError(PS_ERR_IO, false, "failed to read fixts header from skycell file: %s", fileName); + streaksExit("", PS_EXIT_DATA_ERROR); + } + // set up astrometry for conversion from skycell image to sky + pmAstromWCS *wcs = pmAstromWCSfromHeader(header); + if (!wcs) { + psError(PS_ERR_IO, false, "failed to read astrometry from skycell file: %s", fileName); + streaksExit("", PS_EXIT_DATA_ERROR); + } + int naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1"); + int naxis2 = psMetadataLookupS32(NULL, header, "NAXIS2"); + + psFree(header); + + /* now set up our wrapper to the chip astrometry to apply to the whole chip */ + sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, NULL, + sf->diffedPixels->numCols, sf->diffedPixels->numRows); + + // convert corners of skycell to sky coordinates + // compute overlap with the chip using astrometry + psPlane pt[4]; + pt[0].x = 0; + pt[0].y = 0; + + pt[1].x = naxis1 - 1; + pt[1].y = 0; + + pt[2].x = naxis1 - 1; + pt[2].y = naxis2 - 1; + + pt[3].x = 0; + pt[3].y = naxis2 - 1; + + for (int i = 0; i < 4 ; i++ ) { + psSphere sky; + + // convert corner of skycell to sky coordinates + if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) { + psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName); + streaksExit("", PS_EXIT_DATA_ERROR); + } + strkPt p; + // convert to chip coordinates + if (!skyToCell(&p, sf->astrom, sky.r, sky.d)) { + psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName); + streaksExit("", PS_EXIT_DATA_ERROR); + } + pt[i].x = p.x; + pt[i].y = p.y; + } + + psFree(wcs); + + // put the corners in the desired order (see comments below) + + nameCorners(pt); + psString type; + if (pt[0].y < pt[2].y ) { + type = "1"; + } else { + type = "2"; + } +#ifdef DEBUG_PRINT + printf("\nSKYCELL %s Type: %s\n", fileName, type); + for (int i=0; i<4; i++) { + printf("%f %f\n", pt[i].x, pt[i].y); + } +#endif + // Now set the touched pixels + int ymin = fmax(0, pt[1].y ); + int ymax = fmin(pt[3].y + 0.5, sf->diffedPixels->numRows - 1); +#if (DEBUG_PRINT > 1) + printf("\nymin: %d ymax: %d\n", ymin, ymax); +#endif + for (int y = ymin ; y <= ymax; y++) { + int xleft = xLeft(pt, y); + int xright = xRight(pt, y); + + if (xleft < 0) { + xleft = 0; + } + if (xleft > sf->diffedPixels->numCols) { + continue; + } + if (xright < 0) { + continue; + } + if (xright >= sf->diffedPixels->numCols) { + xright = sf->diffedPixels->numCols - 1; + } +#if (DEBUG_PRINT > 1) + printf(" y: %d xleft: %d xright: %d\n", y, xleft, xright); +#endif + + psU8 *rowPixels = sf->diffedPixels->data.U8[y]; + + int n = xright - xleft + 1; + if (n > 0) { + memset(&rowPixels[xleft], 255, n); + } + } +} + +// x as a function of y for the line between two points +// Note: the caller guarentees that the y's of the two points are different +static double xOfY(psPlane *pI, psPlane *pJ, int y) +{ + double m_inv = (double) (pJ->x - pI->x) / (double) (pJ->y - pI->y); + + double x = ((double) y - pI->y) * m_inv + (double) pI->x; + + return x; +} + +static int xLeft(psPlane *pt, int y) +{ + double x_d; + if (y < pt[0].y) { + // left boundary is the line from pt 0 to pt 1 + x_d = xOfY(&pt[0], &pt[1], y); + } else { + // left boundary is the line from pt 0 to pt 3 + x_d = xOfY(&pt[0], &pt[3], y); + } + return (int) (x_d + 0.5); +} + +static int xRight(psPlane *pt, int y) +{ + double x_d; + if (y < pt[2].y) { + // right boundary is the line from pt 1 to pt 2 + x_d = xOfY(&pt[1], &pt[2], y); + + } else { + // right boundary is the line from pt 2 to pt 3 + x_d = xOfY(&pt[2], &pt[3], y); + } + return (int) x_d; +} + +/* + * To compute the overlap of a quadrilateral with a set of + * points we transform the 4 corners to image coordinates + * and name the 4 points of the quad + pt 0 is left most (bottom most corner) + pt 1 is bottom most (right most corner) + pt 2 is right most (top most corner) + pt 3 is top most (left most corner) + +* Then we can divide the quad into 3 regions A B and C based +* on the y coordinate + +* (in the degenerate case of a rectangle aligned to the axes +* we only have 1 region) +* In each of the three regions A, B, and C the test for whether a +* point is contained in the quad on scanline y is simply + + x >= xleft(y) x < xright(y) + +* where xleft(y) and xright(y) are the bounding lines at the given +* y coordinate. + +The following diagrams show some examples. The 4 points are labeled 0 to 3 +The three regions are labeld A B and C, the boundaries between the regions +ocur at pt0.y and pt2.y + +Example 1 + + 3 + C + ---------------2 left boundary: line 0_1 y < pt0.y + B line 0_3 y >= pt0.y + 0-------------- + A right boundary: line 1_2 y < pt2.y + 1 line 2_3 y >= pt.y +************************** +Example 2 + 3 + C + 0---------------- + B + ----------------2 + A + 1 + +************************** +Example 3 + + 3 2 + + B + + 0 1 + +(region A and C are empty in the case where the points form a rectangle) +*/ + +/* + Name the corners + The following algorithm works for points that form a quadrilateral + I think it also works for the situation where 3 points are co-linear + and we have a triangle, but that isn't important for our purposes + +*/ + +static void +nameCorners(psPlane pt[4]) +{ + psPlane pt0, pt1, pt2, pt3; + + /* find pt0 the left most (bottom most) corner */ + int imin = 0; + int min = (int) pt[0].x; + int i; + for (i=1; i<4; i++) { + int x = pt[i].x; + if ((x < min) || + ((x == min) && (pt[i].y < pt[imin].y))) { + imin = i; + min = x; + } + } + // remve pt0 from the list + pt0 = pt[imin]; + for (i=imin; i<3; i++) { + pt[i] = pt[i+1]; + } + + // find the bottom most (right most) corner from the remaining 3 + imin = 0; + min = pt[0].y; + for (i=1; i<3; i++) { + int y = pt[i].y; + if ((y < min) || + ((y == min) && (pt[i].x > pt[imin].x)) ) { + imin = i; + min = y; + } + } + // remve pt1 from the list + pt1 = pt[imin]; + for (i=imin; i<3; i++) { + pt[i] = pt[i+1]; + } + + // now find the right most (top most) of the remaining 2 points + if ((pt[0].x > pt[1].x) || + ((pt[0].x == pt[1].x) && (pt[0].y > pt[1].y))) { + pt2 = pt[0]; + pt3 = pt[1]; + } else { + pt2 = pt[1]; + pt3 = pt[0]; + } + /* write the outputs */ + pt[0] = pt0; + pt[1] = pt1; + pt[2] = pt2; + pt[3] = pt3; +} Index: trunk/magic/remove/src/streaksVersion.c =================================================================== --- trunk/magic/remove/src/streaksVersion.c (revision 26438) +++ trunk/magic/remove/src/streaksVersion.c (revision 26438) @@ -0,0 +1,113 @@ +/** @file streaksVersion.c + * + * @brief + * + * @ingroup streaks + * + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "streaksremove.h" +#include "streaksVersionDefinitions.h" + +#ifndef STREAKS_VERSION +#error "STREAKS_VERSION is not set" +#endif +#ifndef STREAKS_BRANCH +#error "STREAKS_BRANCH is not set" +#endif +#ifndef STREAKS_SOURCE +#error "STREAKS_SOURCE is not set" +#endif + +psString streaksVersion(void) +{ + char *value = NULL; + psStringAppend(&value, "%s@%s", STREAKS_BRANCH, STREAKS_VERSION); + return value; +} + +psString streaksSource(void) +{ + return psStringCopy(STREAKS_SOURCE); +} + +psString streaksVersionLong(void) +{ + psString version = streaksVersion(); // Version, to return + psString source = streaksSource(); // Source + + psStringPrepend(&version, "streaks "); + psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__); + psFree(source); + +#ifdef __OPTIMIZE__ + psStringAppend(&version, " optimised"); +#else + psStringAppend(&version, " unoptimised"); +#endif + + return version; +}; + +bool streaksVersionHeader(psMetadata *header) +{ + PS_ASSERT_METADATA_NON_NULL(header, false); + + psString version = streaksVersion(); // Software version + psString source = streaksSource(); // Software source + + psStringPrepend(&version, "streaks version: "); + psStringPrepend(&source, "streaks source: "); + + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source); + + psFree(version); + psFree(source); + + return true; +} + +bool streaksVersionHeaderFull(psMetadata *header) +{ + PS_ASSERT_METADATA_NON_NULL(header, false); + + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psString history = NULL; // History string + psStringAppend(&history, "streaks at %s", timeString); + psFree(timeString); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history); + psFree(history); + + psLibVersionHeader(header); + psModulesVersionHeader(header); + streaksVersionHeader(header); + + return true; +} + +void streaksVersionPrint(void) +{ + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psLogMsg("streaks", PS_LOG_INFO, "streaks at %s", timeString); + psFree(timeString); + + psString pslib = psLibVersionLong();// psLib version + psString psmodules = psModulesVersionLong(); // psModules version + psString streaks = streaksVersionLong(); // streaks version + + psLogMsg("streaks", PS_LOG_INFO, "%s", pslib); + psLogMsg("streaks", PS_LOG_INFO, "%s", psmodules); + psLogMsg("streaks", PS_LOG_INFO, "%s", streaks); + + psFree(pslib); + psFree(psmodules); + psFree(streaks); + + return; +} Index: trunk/magic/remove/src/streaksVersionDefinitions.h.in =================================================================== --- trunk/magic/remove/src/streaksVersionDefinitions.h.in (revision 26438) +++ trunk/magic/remove/src/streaksVersionDefinitions.h.in (revision 26438) @@ -0,0 +1,8 @@ +#ifndef STREAKS_VERSION_DEFINITIONS_H +#define STREAKS_VERSION_DEFINITIONS_H + +#define STREAKS_VERSION @STREAKS_VERSION@ // SVN version +#define STREAKS_BRANCH @STREAKS_BRANCH@ // SVN branch +#define STREAKS_SOURCE @STREAKS_SOURCE@ // SVN source + +#endif Index: trunk/magic/remove/src/streaksio.c =================================================================== --- trunk/magic/remove/src/streaksio.c (revision 26408) +++ trunk/magic/remove/src/streaksio.c (revision 26438) @@ -121,4 +121,14 @@ sf->transparentStreaks = psMetadataLookupF64(&status, config->arguments, "TRANSPARENT_STREAKS"); + sf->stats = psMetadataAlloc(); + psString statsFileName= psMetadataLookupStr(&status, config->arguments, "STATS"); + if (statsFileName) { + sf->statsFile = fopen(statsFileName, "w"); + if (!sf->statsFile) { + psError(PS_ERR_IO, true, "failed to open stats file %s", statsFileName); + streaksExit("", PS_EXIT_CONFIG_ERROR); + } + } + return sf; } @@ -128,5 +138,5 @@ { freeImages(sf); - psFree(sf->warpedPixels); + psFree(sf->diffedPixels); psFree(sf->tiles); psFree(sf->view); @@ -1181,4 +1191,7 @@ } +// Get the mask values that we need. +// If an input mask file is provided get them from there. +// Otherwise get them from the recipe MASKS void strkGetMaskValues(streakFiles *sfiles) Index: trunk/magic/remove/src/streaksremove.c =================================================================== --- trunk/magic/remove/src/streaksremove.c (revision 26408) +++ trunk/magic/remove/src/streaksremove.c (revision 26438) @@ -12,12 +12,16 @@ static pmConfig *parseArguments(int argc, char **argv); static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll); -static void exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue); -static bool warpedPixel(streakFiles *sfiles, int x, int y); +static void exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue); +static bool diffedPixel(streakFiles *sfiles, int x, int y); static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue); static void writeImages(streakFiles *sf, bool exciseImageCube); static void updateAstrometry(streakFiles *sfiles); static void censorSources(streakFiles *sfiles, psImageMaskType maskStreak); -static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonWarpedPixels, psU16 maskStreak); - +static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonDiffedPixels, psU16 maskStreak); + +// Note: For clarity the flow of this program is in main(). +// There is not a lot of error checking is done in main. +// Until the end, where we might be doing Nebulous operations, called functions exit when an error +// is encountered. int main(int argc, char *argv[]) @@ -26,5 +30,5 @@ psLibInit(NULL); - psTimerStart("STREAKSREMOVE"); + psTimerStart("TOTAL_TIME"); pmConfig *config = parseArguments(argc, argv); @@ -46,4 +50,5 @@ streakFiles *sfiles = openFiles(config, true, argv[0]); setupAstrometry(sfiles); + sfiles->stats = psMetadataAlloc(); // Optionally we can set pixels that are masked to NAN since they couldn't have been @@ -57,12 +62,12 @@ bool exciseAll = false; - // --keepnonwarped is a test and debug mode - bool keepNonWarpedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_WARPED"); - - // we need to check for non warped pixels unless we've been asked not to or the stage is diff - // (By definition pixels in diff images were included in a diff) - bool checkNonWarpedPixels = ! (keepNonWarpedPixels || (sfiles->stage == IPP_STAGE_DIFF) ); - - if (checkNonWarpedPixels ) { + // --keepnondiffed is a test and debug mode + bool keepNonDiffedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_DIFFED"); + + // we need to check for non diffed pixels unless we've been asked not to or the stage is diff + // (By definition pixels in diff images were included in the difference images) + bool checkNonDiffedPixels = ! (keepNonDiffedPixels || (sfiles->stage == IPP_STAGE_DIFF) ); + + if (checkNonDiffedPixels ) { // From magic ICD: // In the raw and detrended images, the pixels which were not @@ -72,11 +77,12 @@ // if no skycells are provided sfiles->exciseAll is set to true - psTimerStart("COMPUTE_WARPED_PIXELS"); - if (! computeWarpedPixels(sfiles) ) { + psTimerStart("COMPUTE_DIFFED_PIXELS"); + if (! computeDiffedPixels(sfiles) ) { // we have no choice to excise all pixels exciseAll = true; } - psF64 cwp_t = psTimerClear("COMPUTE_WARPED_PIXELS"); - psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t); + psF64 cdp_t = psTimerClear("COMPUTE_DIFFED_PIXELS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "COMPUTE_UNDIFFED_PIXELS", PS_META_REPLACE, "time to compute non-diffedpixels", cdp_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to compute diffed pixels: %f\n", cdp_t); } @@ -97,4 +103,10 @@ long totalStreakPixels = 0; + // accumulators for the various timers + psF64 gsp_t = 0; + psF64 enw_t = 0; + psF64 rms_t = 0; + psF64 cs_t = 0; + psF64 wi_t = 0; // Iterate through each component of the input (except for raw images there is only one) do { @@ -126,4 +138,6 @@ StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols, sfiles->inImage->numRows); + gsp_t += psTimerClear("GET_STREAK_PIXELS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "GET_STREAK_PIXELS", PS_META_REPLACE, "", gsp_t); psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS")); @@ -131,19 +145,23 @@ // otherwise it contained an image cube (video cell) which is handled in the if block if (sfiles->inImage->image) { - if (checkNonWarpedPixels) { - psTimerStart("EXCISE_NON_WARPED"); - - // set non-warped pixels and variance to NAN, mask to maskStreak (since the pixel + if (checkNonDiffedPixels) { + psTimerStart("EXCISE_NON_DIFFED"); + + // set non-diffed pixels and variance to NAN, mask to maskStreak (since the pixel // is excised as part of the destreaking process) - exciseNonWarpedPixels(sfiles, sfiles->maskStreak); - - psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED")); + exciseNonDiffedPixels(sfiles, sfiles->maskStreak); + + enw_t += psTimerClear("EXCISE_NON_DIFFED"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "EXCISE_NON_DIFFED", PS_META_REPLACE, "", enw_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non diffed pixels: %f\n", enw_t); } psTimerStart("REMOVE_STREAKS"); - totalStreakPixels += censorPixels(sfiles, pixels, checkNonWarpedPixels, sfiles->maskStreak); - - psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS")); + totalStreakPixels += censorPixels(sfiles, pixels, checkNonDiffedPixels, sfiles->maskStreak); + + rms_t += psTimerClear("REMOVE_STREAKS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "REMOVE_STREAKS", PS_META_REPLACE, "", enw_t); + psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", rms_t); if (nanForRelease) { @@ -167,9 +185,15 @@ } + psTimerStart("CENSOR_SOURCES"); censorSources(sfiles, sfiles->maskStreak); + cs_t += psTimerClear("CENSOR_SOURCES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "CENSOR_SOURCES", PS_META_REPLACE, "", cs_t); + // write the destreaked "temporary" images and the recovery images + psTimerStart("WRITE_IMAGES"); writeImages(sfiles, exciseImageCube); - + wi_t += psTimerClear("WRITE_IMAGES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "WRITE_IMAGES", PS_META_REPLACE, "", wi_t); psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT")); @@ -179,5 +203,11 @@ psFree(streaks); - psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels); + if (exciseAll) { + totalStreakPixels = totalPixels; + } + + psF64 maskedFraction = totalStreakPixels / totalPixels; + psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, maskedFraction * 100); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "MASKED_FRACTION", PS_META_REPLACE, "", maskedFraction); // check the weight and mask files for extra extensions that might be in files @@ -188,10 +218,11 @@ psTimerStart("CLOSE_IMAGES"); - closeImages(sfiles); - - psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES")); - - + psF64 ci_t = psTimerClear("CLOSE_IMAGES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "CLOSE_IMAGES", PS_META_REPLACE, "", ci_t); + + psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", ci_t); + + psTimerStart("REPLICATE_OUTPUTS"); if (!replicateOutputs(sfiles)) { psErrorStackPrint(stderr, "failed to replicate output files"); @@ -199,4 +230,6 @@ exit(PS_EXIT_UNKNOWN_ERROR); } + psF64 ro_t = psTimerClear("REPLICATE_OUTPUTS"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "REPLICATE_OUTPUTS", PS_META_REPLACE, "", ro_t); // NOTE: from here on we can't just quit if something goes wrong. @@ -208,9 +241,9 @@ // swap the instances for the input and output // Note this is a nebulous database operation. No file I/O is performed + psTimerStart("SWAP_INSTANCES"); if (!swapOutputsToInputs(sfiles)) { - // XXX: Now what? // It is up to the program that reverts failed destreak runs to insure that - // any input files that have been swapped are restored and that the de-streaked - // versions are deleted + // any original non-destreaked input files that have been swapped are restored and that the de-streaked + // versions are deleted. psErrorStackPrint(stderr, "failed to swap files"); @@ -218,4 +251,24 @@ // XXX: pick a specific error code for this failure exit(PS_EXIT_UNKNOWN_ERROR); + } + psF64 si_t = psTimerClear("SWAP_INSTANCES"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "SWAP_INSTANCES", PS_META_REPLACE, "", si_t); + } + + psF64 total_time = psTimerClear("TOTAL_TIME"); + psMetadataAddF64(sfiles->stats, PS_LIST_TAIL, "TOTAL_TIME", PS_META_REPLACE, "", total_time); + psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", total_time); + + if (sfiles->statsFile) { + const char *statsMDC = psMetadataConfigFormat(sfiles->stats); + if (!statsMDC || strlen(statsMDC) == 0) { + psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n"); + } else { + fprintf(sfiles->statsFile, "%s", statsMDC); + psFree((void *)statsMDC); + fclose(sfiles->statsFile); + sfiles->statsFile = NULL; + psFree(sfiles->stats); + sfiles->stats = NULL; } } @@ -228,5 +281,4 @@ streaksNebulousCleanup(); pmConfigDone(); - psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE")); psLibFinalize(); @@ -237,5 +289,5 @@ static long -censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonWarpedPixels, psU16 maskStreak) +censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak) { long streakPixels = 0; @@ -245,5 +297,5 @@ if (psImageGet(pixels, x, y)) { ++streakPixels; - if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) { + if (!checkNonDiffedPixels || diffedPixel(sfiles, x, y)) { excisePixel(sfiles, x, y, true, maskStreak); @@ -251,5 +303,5 @@ } else { // This pixel was not included in any warp and has thus already excised - // by exciseNonWarpedPixels + // by exciseNonDiffedPixels } } @@ -268,5 +320,5 @@ fprintf(stderr, "\t-weight WEIGHT.fits: weight file to de-streak\n"); fprintf(stderr, "\t-replace: replace the input images with the output\n"); - fprintf(stderr, "\t-keepnonwarped: do not exise pixels that were not part of difference processing\n"); + fprintf(stderr, "\t-keepnondiffed: do not exise pixels that were not part of difference processing\n"); fprintf(stderr, "\t-transparent val: instead of setting excicsed pixel to NAN add val\n"); fprintf(stderr, "\t-chip_mask MASK.fits: name of mask for chip stage (camera stage mask is passed tih -mask)\n"); @@ -359,8 +411,8 @@ } - if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) { - psArgumentRemove(argnum, &argc, argv); - psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_WARPED", 0, - "skip excising of non warped pixels", true); + if ((argnum = psArgumentGet(argc, argv, "-keepnondiffed"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_DIFFED", 0, + "skip excising of non diffed pixels", true); } @@ -377,5 +429,5 @@ } - // if skycells are not provided then we have to execise all pixels unless -keepnonwarped + // if skycells are not provided then we have to execise all pixels unless -keepnondiffed pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist"); @@ -440,4 +492,11 @@ psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SOURCES", 0, "name of input sources file", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } + + if ((argnum = psArgumentGet(argc, argv, "-stats"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STATS", 0, + "name of input stats file", argv[argnum]); psArgumentRemove(argnum, &argc, argv); } @@ -589,4 +648,5 @@ strkGetMaskValues(sf); + // add the STREAK bit to the mask image pixels setStreakBits(sf->inMask->image, sf->maskStreak); } @@ -721,11 +781,10 @@ sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x]; } - sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = - sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] | newMaskValue; + sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= newMaskValue; } } static void -exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue) +exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue) { int cell_x0 = sfiles->astrom->cell_x0; @@ -749,5 +808,5 @@ } - psU8 *pixels = sfiles->warpedPixels->data.U8[yChip]; + psU8 *pixels = sfiles->diffedPixels->data.U8[yChip]; if (xParity == 1) { @@ -772,10 +831,10 @@ static bool -warpedPixel(streakFiles *sfiles, int x, int y) +diffedPixel(streakFiles *sfiles, int x, int y) { PixelPos chipCoord; if (!CHIP_LEVEL_INPUT(sfiles->stage)) { - // if we're here on a skycell image by definition this pixel was warped + // if we're here on a skycell image by definition this pixel was diffed return true; } @@ -790,12 +849,12 @@ cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y); - if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) { + if (chipCoord.x < 0 || chipCoord.x >= sfiles->diffedPixels->numCols) { return false; } - if (chipCoord.y < 0 || chipCoord.y >= sfiles->warpedPixels->numRows) { + if (chipCoord.y < 0 || chipCoord.y >= sfiles->diffedPixels->numRows) { return false; } - return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false; + return psImageGet(sfiles->diffedPixels, chipCoord.x, chipCoord.y) ? true : false; } Index: trunk/magic/remove/src/streaksremove.h =================================================================== --- trunk/magic/remove/src/streaksremove.h (revision 26408) +++ trunk/magic/remove/src/streaksremove.h (revision 26438) @@ -64,4 +64,6 @@ sFile *inSources; sFile *outSources; + psMetadata *stats; + FILE *statsFile; psString class_id; pmFPAfile *inAstrom; @@ -71,5 +73,5 @@ pmChip *chip; // current chip pmCell *cell; // current cell - psImage *warpedPixels; + psImage *diffedPixels; psVector *tiles; double transparentStreaks; @@ -88,5 +90,5 @@ extern void cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell); -extern bool computeWarpedPixels(streakFiles *sf); +extern bool computeDiffedPixels(streakFiles *sf); extern void streaksExit(psString, int); extern ippStage parseStage(psString); Index: trunk/magic/remove/src/warpedpixels.c =================================================================== --- trunk/magic/remove/src/warpedpixels.c (revision 26408) +++ (revision ) @@ -1,330 +1,0 @@ -#include "streaksremove.h" -#include "pmAstrometryWCS.h" - -// #define DEBUG_PRINT 1 - -static void addTouchedPixels(streakFiles *sf, psString fileName); -static void nameCorners(psPlane pt[4]); -static int xLeft(psPlane *pt, int y); -static int xRight(psPlane *pt, int y); - -bool -computeWarpedPixels(streakFiles *sf) -{ - bool status; - pmConfig *config = sf->config; - - psArray *skycells = psMetadataLookupPtr(&status, config->arguments, "SKYCELLS"); - if (skycells == NULL) { - return false; - } - - int n = psArrayLength(skycells); - - psRegion *bounds = pmChipPixels(sf->chip); - - int width = bounds->x1 - bounds->x0; - int height = bounds->y1 - bounds->y0; - - psFree(bounds); - - sf->warpedPixels = psImageAlloc(width, height, PS_TYPE_U8); - psImageInit(sf->warpedPixels, 0); - - for (int i=0; iinAstrom->fpa->toSky->type != PS_PROJ_DIS) { - pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001); - } - psFits *fits = psFitsOpen("warpedpixels.fits", "w"); - - psFitsWriteImage(fits, header, sf->warpedPixels, 0, NULL); - psFitsClose(fits); - psFree(header); - } - return true; -} - -static void -addTouchedPixels(streakFiles *sf, psString fileName) -{ - // read fits header - psFits *fits = psFitsOpen(fileName, "r"); - if (!fits) { - psError(PS_ERR_IO, false, "failed to open skycell file: %s", fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - psMetadata *header = psFitsReadHeader(NULL, fits); - psFitsClose(fits); - if (!header) { - psError(PS_ERR_IO, false, "failed to read fixts header from skycell file: %s", fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - // set up astrometry for conversion from skycell image to sky - pmAstromWCS *wcs = pmAstromWCSfromHeader(header); - if (!wcs) { - psError(PS_ERR_IO, false, "failed to read astrometry from skycell file: %s", fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - int naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1"); - int naxis2 = psMetadataLookupS32(NULL, header, "NAXIS2"); - - psFree(header); - - /* now set up our wrapper to the chip astrometry to apply to the whole chip */ - 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 - // compute overlap with the chip using astrometry - psPlane pt[4]; - pt[0].x = 0; - pt[0].y = 0; - - pt[1].x = naxis1 - 1; - pt[1].y = 0; - - pt[2].x = naxis1 - 1; - pt[2].y = naxis2 - 1; - - pt[3].x = 0; - pt[3].y = naxis2 - 1; - - for (int i = 0; i < 4 ; i++ ) { - psSphere sky; - - // convert corner of skycell to sky coordinates - if (!pmAstromWCStoSky(&sky, wcs, &pt[i])) { - psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - strkPt p; - // convert to chip coordinates - if (!skyToCell(&p, sf->astrom, sky.r, sky.d)) { - psErrorStackPrint(stderr, "failed to convert pt %d of %s to sky coords", i, fileName); - streaksExit("", PS_EXIT_DATA_ERROR); - } - pt[i].x = p.x; - pt[i].y = p.y; - } - - psFree(wcs); - - // put the corners in the desired order (see comments below) - - nameCorners(pt); - psString type; - if (pt[0].y < pt[2].y ) { - type = "1"; - } else { - type = "2"; - } -#ifdef DEBUG_PRINT - printf("\nSKYCELL %s Type: %s\n", fileName, type); - for (int i=0; i<4; i++) { - printf("%f %f\n", pt[i].x, pt[i].y); - } -#endif - // Now set the touched pixels - int ymin = fmax(0, pt[1].y ); - int ymax = fmin(pt[3].y + 0.5, sf->warpedPixels->numRows - 1); -#if (DEBUG_PRINT > 1) - printf("\nymin: %d ymax: %d\n", ymin, ymax); -#endif - for (int y = ymin ; y <= ymax; y++) { - int xleft = xLeft(pt, y); - int xright = xRight(pt, y); - - if (xleft < 0) { - xleft = 0; - } - if (xleft > sf->warpedPixels->numCols) { - continue; - } - if (xright < 0) { - continue; - } - if (xright >= sf->warpedPixels->numCols) { - xright = sf->warpedPixels->numCols - 1; - } -#if (DEBUG_PRINT > 1) - printf(" y: %d xleft: %d xright: %d\n", y, xleft, xright); -#endif - - psU8 *rowPixels = sf->warpedPixels->data.U8[y]; - - int n = xright - xleft + 1; - if (n > 0) { - memset(&rowPixels[xleft], 255, n); - } - } -} - -// x as a function of y for the line between two points -// Note: the caller guarentees that the y's of the two points are different -static double xOfY(psPlane *pI, psPlane *pJ, int y) -{ - double m_inv = (double) (pJ->x - pI->x) / (double) (pJ->y - pI->y); - - double x = ((double) y - pI->y) * m_inv + (double) pI->x; - - return x; -} - -static int xLeft(psPlane *pt, int y) -{ - double x_d; - if (y < pt[0].y) { - // left boundary is the line from pt 0 to pt 1 - x_d = xOfY(&pt[0], &pt[1], y); - } else { - // left boundary is the line from pt 0 to pt 3 - x_d = xOfY(&pt[0], &pt[3], y); - } - return (int) (x_d + 0.5); -} - -static int xRight(psPlane *pt, int y) -{ - double x_d; - if (y < pt[2].y) { - // right boundary is the line from pt 1 to pt 2 - x_d = xOfY(&pt[1], &pt[2], y); - - } else { - // right boundary is the line from pt 2 to pt 3 - x_d = xOfY(&pt[2], &pt[3], y); - } - return (int) x_d; -} - -/* - * To compute the overlap of a quadrilateral with a set of - * points we transform the 4 corners to image coordinates - * and name the 4 points of the quad - pt 0 is left most (bottom most corner) - pt 1 is bottom most (right most corner) - pt 2 is right most (top most corner) - pt 3 is top most (left most corner) - -* Then we can divide the quad into 3 regions A B and C based -* on the y coordinate - -* (in the degenerate case of a rectangle aligned to the axes -* we only have 1 region) -* In each of the three regions A, B, and C the test for whether a -* point is contained in the quad on scanline y is simply - - x >= xleft(y) x < xright(y) - -* where xleft(y) and xright(y) are the bounding lines at the given -* y coordinate. - -The following diagrams show some examples. The 4 points are labeled 0 to 3 -The three regions are labeld A B and C, the boundaries between the regions -ocur at pt0.y and pt2.y - -Example 1 - - 3 - C - ---------------2 left boundary: line 0_1 y < pt0.y - B line 0_3 y >= pt0.y - 0-------------- - A right boundary: line 1_2 y < pt2.y - 1 line 2_3 y >= pt.y -************************** -Example 2 - 3 - C - 0---------------- - B - ----------------2 - A - 1 - -************************** -Example 3 - - 3 2 - - B - - 0 1 - -(region A and C are empty in the case where the points form a rectangle) -*/ - -/* - Name the corners - The following algorithm works for points that form a quadrilateral - I think it also works for the situation where 3 points are co-linear - and we have a triangle, but that isn't important for our purposes - -*/ - -static void -nameCorners(psPlane pt[4]) -{ - psPlane pt0, pt1, pt2, pt3; - - /* find pt0 the left most (bottom most) corner */ - int imin = 0; - int min = (int) pt[0].x; - int i; - for (i=1; i<4; i++) { - int x = pt[i].x; - if ((x < min) || - ((x == min) && (pt[i].y < pt[imin].y))) { - imin = i; - min = x; - } - } - // remve pt0 from the list - pt0 = pt[imin]; - for (i=imin; i<3; i++) { - pt[i] = pt[i+1]; - } - - // find the bottom most (right most) corner from the remaining 3 - imin = 0; - min = pt[0].y; - for (i=1; i<3; i++) { - int y = pt[i].y; - if ((y < min) || - ((y == min) && (pt[i].x > pt[imin].x)) ) { - imin = i; - min = y; - } - } - // remve pt1 from the list - pt1 = pt[imin]; - for (i=imin; i<3; i++) { - pt[i] = pt[i+1]; - } - - // now find the right most (top most) of the remaining 2 points - if ((pt[0].x > pt[1].x) || - ((pt[0].x == pt[1].x) && (pt[0].y > pt[1].y))) { - pt2 = pt[0]; - pt3 = pt[1]; - } else { - pt2 = pt[1]; - pt3 = pt[0]; - } - /* write the outputs */ - pt[0] = pt0; - pt[1] = pt1; - pt[2] = pt2; - pt[3] = pt3; -}