Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 25404) @@ -25,4 +25,6 @@ err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell); +// fprintf(stderr, "Got through burn_check\n"); + /* Expand the masks; burns asymmetrically, stars symmetrically */ // err = grow_mask(nx, ny, NX, mbuf, YMASK_GROW, YMASK_GROW, cell->nstar, cell->star); @@ -30,6 +32,8 @@ err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW, MASK_SAT_HALO, cell->nburn, cell->burn); +// fprintf(stderr, "Got through grow mask for burns\n"); err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW, MASK_STAR_HALO, cell->nstar, cell->star); +// fprintf(stderr, "Got through grow mask for stars\n"); /* Fix up all the burns */ @@ -37,7 +41,9 @@ if(!cell->burn[k].burned) continue; /* Fit the trail */ +// fprintf(stderr, "Fitting trail %d\n", k); fit_trail(nx, ny, NX, imbuf, mbuf, cell->burn+k, 1, cell->sky+cell->bias, cell->rms, BURN_PWR); /* Subtract the fit from the image */ +// fprintf(stderr, "Subtracting trail %d\n", k); if(!cell->burn[k].fiterr) sub_fit(nx, ny, NX, buf, cell->burn+k, 1); } @@ -218,7 +224,102 @@ /****************************************************************/ /* burn_test() tests whether a box is burned or not */ +/* It's trying to balance the flux of above-below to see whether + * there's a net burn above. It has to cope with the case that the box is + * so low there is no above in which case it compares center-side. + * It's somewhat rudimentary as indicated by the FIXME's, and the original + * version had some real logical problems. + */ STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms, MTYPE *mask, OBJBOX *box) { + int i, j, nburn, nref, nmed; + int y0, y1, ymid, xburn, xref, dx=0, dy=0; + +/* FIXME: needs a test for flat-toppedness as well as trail... */ +/* FIXME: needs a test for pure, massive saturation */ + +/* Center line */ + ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4; + +/* Irrelevant, too near the top to detect a burn anyway */ + if(box->ey > ny-FIT_EDGE-3) { + box->diff = -1; + box->burned = -1; +/* But looks like a bad one which will leave behind persistence */ + if(box->ey-box->sy+1 > 4) box->burned = 1; + return(0); + +/* Too near the bottom, do a side-side test */ + } else if(box->sy < FIT_EDGE+3) { + dy = 0; + y0 = ymid + MAX(box->ey-ymid, FIT_EDGE); + y1 = MIN(y0+100, ny-1); + + if(2*box->ex - box->sx + 1 < nx-1) { /* Work right? */ + dx = box->ex - box->sx + 1; + + } else if(2*box->sx - box->ex - 1 > 0) { /* Work left? */ + dx = -(box->ex - box->sx + 1); + + } else { /* No room! Oh no! */ + box->diff = -1; + box->burned = -1; +/* But looks like a bad one which will leave behind persistence */ + if(box->ey-box->sy+1 > 4) box->burned = 1; + return(0); + } + +/* Do a top-bottom test for excess flux */ + } else { + dx = 0; + dy = MAX(ymid-box->sy, box->ey-ymid); + if(dy < FIT_EDGE) dy = FIT_EDGE; + y0 = ymid + dy; + y1 = y0 + MIN(ny-1-y0, 2*ymid-y0-1); + if(y1-y0 > 100) y1 = y0 + 100; + } + +/* Do the test; get the median of the averages across x */ + nmed = 0; + for(j=y0; jsx; i<=box->ex; i++) { + if(mask[i+j*NX] == MASK_NONE) { + xburn += data[i+j*NX]; + nburn++; + } + if(dx != 0) { /* Side to side */ + if(mask[i+dx+j*NX] == MASK_NONE) { + xref += data[i+dx+j*NX]; + nref++; + } + } else { /* top-bottom */ + if(mask[i+(2*ymid-j)*NX] == MASK_NONE) { + xref += data[i+(2*ymid-j)*NX]; + nref++; + } + } + } + if(nburn > 0) xburn /= nburn; + if(nref > 0) xref /= nref; + if(nref > 0 && nburn > 0) median_buf[nmed++] = xburn - xref; + } + box->diff = int_median(nmed, median_buf); + box->burned = (box->diff > 0.3*rms); + return(0); +} + +#ifdef ORIG_BURN_TEST +// This has a lot of problems -- chief among which is what happens with +// a huge kite-shaped bad region. It's fundamentally trying to balance +// the flux of above-below to see whether there's a net burn above. But +// it has to cope with the case that the box is so low there is no above +// in which case it compares center-side. +/****************************************************************/ +/* burn_test() tests whether a box is burned or not */ +STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms, + MTYPE *mask, OBJBOX *box) +{ int i, j, n, nrow, nmed; int y0, y1, ymid, xsum, dx=0; @@ -229,4 +330,5 @@ /* Center line */ ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4; + /* Starting point offset */ y0 = MAX(ymid-box->sy, box->ey-ymid); @@ -259,4 +361,5 @@ nrow = 0; for(i=box->sx; i<=box->ex; i++) { +// IS THIS THE PROBLEM? xsum += data[i+(ymid+y0+j)*NX]*(mask[i+(ymid+y0+j)*NX] == MASK_NONE); if(y1 < ny) { @@ -279,2 +382,3 @@ return(0); } +#endif Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnparams.h =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnparams.h (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnparams.h (revision 25404) @@ -55,3 +55,5 @@ EXTERN int EXPIRE_TRAIL_TIME; /* Expire trails after this interval */ +EXTERN int PERSIST_RETAIN; /* Retain bad-slope persistence fits */ + #endif /* _INCLUDED_burnparams_ */ Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 25404) @@ -100,4 +100,6 @@ EXPIRE_TRAIL_TIME = 2000; /* Expire a persist after this [sec] */ + PERSIST_RETAIN = 0; /* Retain persists with bad slopes? */ + /* Parse the args */ cellxy = -1; @@ -175,4 +177,8 @@ persistfitsfile = argv[i] + 7; +/* Keep persistence streaks which had a bad slope? */ + } else if(strncmp(argv[i], "persist=", 8) == 0) {/* persist={t|f} */ + PERSIST_RETAIN = argv[i][8] == 'y' || argv[i][8] == '1' || argv[i][8] == 't'; + /* Output file for PSF gallery */ } else if(strncmp(argv[i], "psf=", 4) == 0) { /* psf=fname */ @@ -286,4 +292,19 @@ CONCAT_FITS = argv[i][10] == 'n' || argv[i][10] == '0' || argv[i][10] == 'f'; + +/* Set the e/ADU for the noise gate that cells must pass */ + } else if(strncmp(argv[i], "EADU=", 5) == 0) { /* EADU=e */ + if(sscanf(argv[i]+5, "%lf", &MIN_EADU) != 1) { + fprintf(stderr, "\rerror: cannot get e/ADU size from `%s'\n", argv[i]); + exit(EXIT_FAILURE); + } + + + /* Set the read noise for the noise gate that cells must pass */ + } else if(strncmp(argv[i], "RN=", 3) == 0) { /* EADU=e */ + if(sscanf(argv[i]+3, "%d", &MAX_READ_NOISE) != 1) { + fprintf(stderr, "\rerror: cannot get read noise size from `%s'\n", argv[i]); + exit(EXIT_FAILURE); + } /* Quiet? */ Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 25404) @@ -97,4 +97,5 @@ #define FIT_TOP_ERROR 2 /* Saturation extends to top: no points */ #define FIT_SLOPE_ERROR 3 /* Unreasonable fit */ +#define FIT_ALL_GONE 4 /* No column survived as significant */ #define FIT_EXPIRED 9 /* Don't carry any more as persistent */ Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c (revision 25404) @@ -115,4 +115,5 @@ if(VERBOSE & VERB_BOXGROW) { printf("box: %d %d %d %d\n", box->sx, box->sy, box->ex, box->ey); + fflush(stdout); } Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 25404) @@ -53,4 +53,12 @@ power laws; downward burns as exponentials. + The fit to downward, persistence burns may have an unreasonable "slope" + (meaning exponential sign), increasing away from the burn origin. + These are discarded unless "persist=t" is invoked. It is possible + for the fit to be fooled by uncataloged stars, so the persistence + can be kept in the output file with "persist=t" (although no fit + correction is applied) and they will propagate from input to output + until they finally achieve a legal fit with negligible amplitude. + Burntool also identifies really blasted areas which are saturated from top to bottom. These cannot be fitted, but are carried along for @@ -231,4 +239,8 @@ 'in' takes precedence. + persist={t|f} + Retain persistence streaks whose fit slope turned out to be + unreasonable. + out=fname Output file for burn streaks @@ -290,4 +302,15 @@ expire=N Retire a blasted burn after N seconds + + EADU=x + Set the minimum e/ADU for acceptable cell sky noise + (default 0.3). Should not normally be changed! + RN=x + Set the maximum read noise for acceptable cell sky noise + (default 20). Should not normally be changed! + + Any cell whose noise variance in the sky is greater than + sky/EADU+RN*RN is rejected out of hand, and skipped for + further processing. quiet={t|f} @@ -332,4 +355,5 @@ BUGS: 090224: Still in development + 090810: Squished a few memory bugs SEE ALSO: Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/persist_fits.c =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/persist_fits.c (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/burntool/persist_fits.c (revision 25404) @@ -222,6 +222,14 @@ for(k=0; k 0) { + boxbuf[nbox].zero = (double *)calloc(boxbuf[nbox].nfit, sizeof(double)); + boxbuf[nbox].xfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int)); + boxbuf[nbox].yfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int)); + for(i=0; i 0) zk = box[k].zero[box[k].nfit/2]; for(jp=kp+1; jp 0) zj = box[j].zero[box[j].nfit/2]; if(ABS(yctr[jp]-yctr[kp]) > DIFFERENT_STREAK) { /* Trim back the feebler streak */ @@ -238,6 +241,15 @@ /* First: patched up persists */ for(k=0; k= FIT_MAX_SLOPE || slope < FIT_MIN_SLOPE) { + box->slope = slope; + box->nfit = 0; box->fiterr = FIT_SLOPE_ERROR; return(-1); @@ -283,7 +285,13 @@ } /* Update fit ranges */ - box->sxfit = box->xfit[0]; - box->exfit = box->xfit[box->nfit-1]; - + if(box->nfit > 0) { + box->sxfit = box->xfit[0]; + box->exfit = box->xfit[box->nfit-1]; + } else { + box->sxfit = xs; + box->exfit = xs - 1; + box->fiterr = FIT_ALL_GONE; + } +// if(box->nfit <= 0) fprintf(stderr, "WHOA, got nfit = 0\n"); return(0); } Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/libfh/fh.c =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/libfh/fh.c (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/libfh/fh.c (revision 25404) @@ -1039,4 +1039,5 @@ newval[strlen(newval) - 1] = '\0'; fh_set_str(hu, idx, name, newval, comment); + free(newval); return; } @@ -1519,4 +1520,29 @@ i++; } + return FH_SUCCESS; +} + +static void +pad_header(HeaderUnit hu, int n) +{ + int i; + double idx = 900000000; + + for (i = 0; i < n; i++) + { + fh_set_card(hu, idx++, FH_RESERVE); + } +} + +fh_result +fh_copy(HeaderUnit hu, const HeaderUnit source_) +{ + fh_result result; + int reserve; + + result = fh_merge(hu, source_); + if (result != FH_SUCCESS) return result; + reserve = fh_get_reserve(source_); + if (reserve) pad_header(hu, reserve); return FH_SUCCESS; } @@ -1757,4 +1783,21 @@ } +int +fh_get_reserve(HeaderUnit hu) +{ + HeaderUnitStruct* list = FH_HU(hu); + + if (!list) + { + log_error("invalid argument to fh_get_reserve()"); + return -1; + } + + if (list->reserve) + return list->reserve; + else + return list->reserve_found; +} + fh_result fh_rewrite(HeaderUnit hu) Index: /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/libfh/fh.h =================================================================== --- /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/libfh/fh.h (revision 25403) +++ /branches/eam_branches/20090715/extsrc/gpcsw/gpcsrc/fits/libfh/fh.h (revision 25404) @@ -174,4 +174,9 @@ * is used to create a new FITS header. NOT for use with fh_rewrite(). */ +int fh_get_reserve(HeaderUnit hu); +/* + * Get the current setting or the number of reserve cards found in + * a header unit. Returns -1 if hu is invalid. + */ fh_result fh_validate(HeaderUnit hu); /* fh_validate.c; for use by fhtool.c */ @@ -386,4 +391,5 @@ double fh_idx(HeaderUnit hu); /* idx of the last card returned by fh_next */ fh_result fh_merge(HeaderUnit hu, const HeaderUnit source); /* source unchanged */ +fh_result fh_copy(HeaderUnit hu, const HeaderUnit source); /* fh_merge+keep reserve */ /* --------------------------------------------------------- @@ -459,4 +465,3 @@ */ - #endif /* _INCLUDED_fh */