Index: trunk/ppImage/src/Makefile.am =================================================================== --- trunk/ppImage/src/Makefile.am (revision 35612) +++ trunk/ppImage/src/Makefile.am (revision 35685) @@ -66,5 +66,6 @@ ppImageAddNoise.c \ ppImageRandomGaussian.c \ - ppImageAuxiliaryMask.c + ppImageAuxiliaryMask.c \ + ppImageBurntoolApply.c CLEANFILES = *~ Index: trunk/ppImage/src/burnparams.h =================================================================== --- trunk/ppImage/src/burnparams.h (revision 35685) +++ trunk/ppImage/src/burnparams.h (revision 35685) @@ -0,0 +1,59 @@ +/* A global buffer into which routines can accumulate detections */ +#ifndef _INCLUDED_burnparams_ +#define _INCLUDED_burnparams_ + +/* '#define EXTERN' in just one file, burntool.c, to declare variables */ +#ifndef EXTERN +#define EXTERN extern +#endif + +#define MAXBURN 10000 +EXTERN OBJBOX boxbuf[MAXBURN]; + +EXTERN int *median_buf; /* Generic buffer for integer medians */ +EXTERN int nmedian_buf; + +EXTERN MTYPE *mbuf; /* Buffer for object mask */ +EXTERN int nmbuf; + +EXTERN MTYPE *msbuf; /* Buffer for vetoed star mask */ +EXTERN int nmsbuf; + +EXTERN DTYPE *imbuf; /* Copy of cell data w/o BZERO */ +EXTERN int nimbuf; + +EXTERN int VERBOSE; /* Verbosity level */ + +EXTERN int BZERO; +EXTERN int USHORT_BIAS; /* Bias level restored to ushort stamps */ + +EXTERN int MAX_READ_NOISE; /* Maximum believable read noise (ADU) */ +EXTERN double MIN_EADU; /* Minimum believable e/ADU */ +EXTERN int SAT4SURE; /* Ignore pixels above for noise estimate */ +EXTERN double MIN_BLAST_PASS; /* Allow blasted cells if they have BIG satfrac */ +EXTERN double MAX_BLAST_PASS; /* But not if it's all wiped out! */ + +EXTERN int BURN_THRESH ; /* Threshold for onset of burning */ +EXTERN int TRAIL_THRESH ; /* Trailing might go this low */ +EXTERN int MAX_THRESH ; /* Possibly trailing stars? */ +EXTERN int STAR_THRESH ; /* Threshold for star above sky */ +EXTERN double STAR_FRAC ; /* Fraction to follow star profile */ +EXTERN int PSF_THRESH ; /* Threshold for a star to be a PSF */ + +//EXTERN double XMASK_GROW ; /* Growth of burned boxes in x dir */ +//EXTERN double YMASK_GROW ; /* Growth of burn/star in y dir */ + +EXTERN double BMASK_GROW ; /* Growth of burned boxes in size */ +EXTERN double RMASK_GROW ; /* Growth of burn/star in diameter */ + +EXTERN int MIN_PSF_SIZE ; /* Min box size for stamp selection */ +EXTERN int PSF_CTR_TOL ; /* Choose max or box ctr for stamp */ +EXTERN int CONCAT_FITS; /* Write concat FITS for PSF? (else 3D) */ +EXTERN int MAX_PSF_PER_CELL ; /* Max number of PSF stars accepted per cell */ + +EXTERN double NEGLIGIBLE_TRAIL; /* Don't sweat less than this * sigma */ +EXTERN int EXPIRE_TRAIL_TIME; /* Expire trails after this interval */ + +EXTERN int PERSIST_RETAIN; /* Retain bad-slope persistence fits */ + +#endif /* _INCLUDED_burnparams_ */ Index: trunk/ppImage/src/burntool.h =================================================================== --- trunk/ppImage/src/burntool.h (revision 35685) +++ trunk/ppImage/src/burntool.h (revision 35685) @@ -0,0 +1,229 @@ +#ifndef _INCLUDED_burntool_ +#define _INCLUDED_burntool_ + +/* Burn correction routines: + ------------------------- + * 1. Identify burned spots from saturation threshold + ? 2. Identify burns and trails morphologically? + ? 2.1 Look for flat-topped stars + ? 2.2 Look for asymmetric, smooth, declining trails? + ? 2.3 Use y overscan? + ? 2.4 Bin in x? + * 3. Mask stars and other cruft? + * 4. Fit up and down row functions + * 5. Read burn time list + * 6. Write updated burn time list + * 8. Undo subtraction + 7. Create and save burn correction FITS table +*/ +/* Algorithm: + ---------- + 1. Find new burned spots + 2. Create new burn entries: x1,x2,y,time=now + 3. For each burn: + if t = now detrail up + else deburn down + if fit negligible delete burn + else save burn fits + */ +/* ToDo: + ----- + * clean up fit start points, esp for persistence + * change burn= and persist= to trailin= and trailout= + * trim edges of fits + * implement MASK_SAT + > test for garbage fits (positive slope, etc) + * test nearly saturated stars for burn. + * undo subtraction (almost perfect -- occasional pixels and baddies) + * star gather into 3D FITS + * test on donuts + * stamp extraction center on box, not max + * Add back nominal bias; output as ushort + * man page... ongoing... + * Scale uses peak/sum? + * Change psf coords to center; save in header not pixels + x better trail function? A tuned ln(x+A) is better, but too hard. + * why various restorations imperfect? edge of cell? + * cell mask argument +... why the honker star not ID'ed as a burn? various other gotchas... + * age away old persistence streaks; merge persistence streaks + * PSF info on median psf + * clean up sky/rms calculation with 5*quartile clipping + * change VERBOSE to a bitmask for access to various portions of the code + * add access to a lot of the tuning variables... + * update man page + tuning, tuning, tuning... + add correction FITS table to MEF +*/ + +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define NINT(x) (x<0?(int)((x)-0.5):(int)((x)+0.5)) +#define ABS(a) (((a) > 0) ? (a) : -(a)) + +// #define STATIC static /* To make function declarations static */ +#define STATIC /* To make function declarations global */ + +#define MAXCELL 64 /* Max cells in an OTA */ +#define MAXSIZE 2048 /* Maximum vertical cell size */ + +#define STAR_RADIUS 4 /* Radius over which a star ctr must be max */ +#define SKY_MARG 3 /* Horiz offset from burn area for sky */ +#define FIT_EDGE 10 /* How far beyond saturation to start fit? */ +#define Y_SCALE 0.01 /* Scale factor for y in fits */ + +typedef signed short int IMTYPE;/* Data type of image */ + +typedef int DTYPE; /* Data type of data copy */ +#define NODATA 0 /* Marker for *No Data */ + +/* Mask codes: keep in order of severity for psf_select()! */ +typedef int MTYPE; /* Data type of mask */ +#define MASK_NONE 0 /* Unmasked pixel (MUST be 0) */ +#define MASK_STAR_HALO 1 /* Extended halo box from star */ +#define MASK_SAT_HALO 2 /* Extended halo box from grow_mask() */ +#define MASK_CTR 3 /* Center mask from star_detect() */ +#define MASK_SAT 4 /* Center mask of saturated pixels */ + +/* Enumeration of function choices */ +#define FUNC_NONE 0 /* No associated function */ +#define BURN_PWR 1 /* Power law */ +#define BURN_EXP 2 /* Exponential */ +#define BURN_BLASTED 3 /* Blasted top to bottom: flag only for IPP */ +#define BURN_POSSLOPE 4 /* Positive slope fit (bad) but significant */ +#define PSF_STAR 9 /* Unfitted: good psf star */ + +/* Fit error codes */ +#define FIT_ERROR 1 /* linearfit failed */ +#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 */ + +/* Fit parameters */ +#define FIT_MIN_SLOPE -5.0 /* minimum slope which is a credible fit */ +#define FIT_MAX_SLOPE 0.0 /* maximum slope which is a credible fit */ + +/* Verbosity bits */ +#define VISTAMARKER "/tmp/markem.pro" /* Write Vista pro to mark stars?*/ +#define VERB_NORM 0x0001 /* Normal, verbose output */ +#define VERB_DETECT 0x0002 /* Dump detection process */ +#define VERB_PSFSEL 0x0004 /* Dump out PSF selection process */ +#define VERB_FIT 0x0008 /* Dump fit progress */ +#define VERB_FITPROF 0x0010 /* Dump fit profiles */ +#define VERB_VISTA 0x0020 /* Write vista markers as /tmp/markem.pro */ +#define VERB_BOXGROW 0x0040 /* Dump box growth diagnostics */ +#define VERB_MASK 0x0080 /* Write mask in place of corrected image */ + +/* Description of a potentially burned trail */ +typedef struct obj_box { + int cell; /* what cell is this one in? */ + int time; /* PON time when it was created */ + int sx; /* left corner */ + int sy; /* bottom corner */ + int ex; /* right corner */ + int ey; /* top corner */ + int cx; /* center x (position of max) */ + int cy; /* center y (position of max) */ + int max; /* max data value (above sky) */ + int y0m; /* min y value at sx */ + int y0p; /* max y value at sx */ + int y1m; /* min y value at ex */ + int y1p; /* max y value at ex */ + int x0m; /* min x value at sy */ + int x0p; /* max x value at sy */ + int x1m; /* min x value at ey */ + int x1p; /* max x value at ey */ + int sat; /* saturated (i.e. bigger than BURNTHRESH)? */ + int func; /* what are we going to do about it? */ + int diff; /* median diff of up minus down (sum for stars) */ + int up; /* does it trail up or down? (stamp sx for stars) */ + int y0; /* y origin for the fit (stamp sy for stars) */ + int burned; /* do we think it's left a burn trail? */ + int fiterr; /* error fitting the trail? */ + int sxfit; /* starting column for fits */ + int exfit; /* ending column for fits */ + int nfit; /* how many columns were corrected? */ + IMTYPE *stamp; /* postage stamp of this object */ + int *xfit; /* x of each value of the start of correction */ + int *yfit; /* y value of the start of correction */ + double slope; /* slope of fit (stamp sum/max) */ + double *zero; /* zero of fit for each of the columns */ +} OBJBOX; + +/* Info for an entire cell */ +typedef struct cell_info { + int cell; /* Cell number */ + int bias; /* Bias level */ + int sky; /* Sky level */ + int rms; /* RMS in the sky */ + int time; /* PON time of this cell */ + double satfrac; /* Fraction of SAT4SURE saturated pixels */ + int nburn; /* Number of trails left by sat stars */ + OBJBOX *burn; /* Stars which we think left a trail */ + int npersist; /* Number of old persistence streaks */ + OBJBOX *persist; /* Persistent streaks */ + int nstar; /* Number of stars */ + OBJBOX *star; /* Stars, harmless we think */ +} CELL; + +/* Prototypes */ +STATIC int mem_init(int nx, int ny, int NX, int NY); +STATIC int burn_fix(int nx, int ny, int stride, int NY, IMTYPE *buf, + CELL *cell, int cellnum); +STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms, + MTYPE *mask, OBJBOX *box); +STATIC int burn_restore(int nx, int ny, int NX, IMTYPE *buf, CELL *cell); +STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell); +STATIC int persist_read(CELL *cell, const char *infile, int apply); +STATIC int persist_write(CELL *cell, const char *outfile); +STATIC int persist_fix(int nx, int ny, int stride, IMTYPE *buf, CELL *cell); +STATIC int persist_merge(CELL *cell); + +//fh_result persist_fits_read(CELL *cell, const char *filename, int apply); +//fh_result persist_fits_write(CELL *cell, HeaderUnit phu); +//fh_result persist_fits_remove_tables(HeaderUnit phu_in, const char *fileout); + +STATIC int star_detect(int nx, int ny, int NX, int NY, DTYPE *data, + MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum); +STATIC int burn_check(int nx, int ny, int stride, int NY, DTYPE *buf, + MTYPE *mask, CELL *cell); + +STATIC int cell_stats(int nx, int ny, int NX, int NY, DTYPE *data, CELL *cell); +STATIC int burn_blab(CELL *cell); +STATIC int persist_blab(CELL *cell); +STATIC int vista_marker(CELL *cell, char *fname); + +STATIC int fit_trail(int nx, int ny, int NX, DTYPE *data, MTYPE *mask, + OBJBOX *box, int up, int sky, int rms, int fitfunc); +STATIC int sub_fit(int nx, int ny, int NX, IMTYPE *buf, OBJBOX *box, int sign); + +STATIC int psf_select(int nx, int ny, int NX, MTYPE *mask, DTYPE *buf, + int nbox, OBJBOX *box, int size, int sky); +STATIC int psf_write(int nx, int ny, CELL *OTA, int otanum, const char *psffile); +STATIC int psf_write_stats(int nx, int ny, CELL *OTA, int otanum, const char *statfile, int psfavg); +STATIC int psf_stats(int nx, int ny, IMTYPE *data, int bias, + double *fwhm, double *q); + +STATIC int grow_box(int nx, int ny, int NX, DTYPE *data, int thr, OBJBOX *box); +STATIC int grow_mask(int nx, int ny, int NX, DTYPE *mask, double xfac, + double yfac, int maskval, int nbox, OBJBOX *box); +STATIC int local_max(int i, int j, int r, int nx, int ny, int NX, DTYPE *data); + +STATIC int wlinearfit(int npt, double *x, double *y, + double *w, double *a, double *b); +STATIC int linearrms(int npt, double *x, double *y, double a, double b, double *rms); +STATIC int int_median(int n, int *key); +STATIC double double_median(int n, double *key); + +STATIC int qsort_int(int n, int *key); +STATIC double qsort_dbl(int n, double *key); + +STATIC void syntax(const char *prog); + +int write_2dfits(int nx, int ny, int sx, int sy, IMTYPE *data, int fd); +int write_3dhdr(int nx, int ny, int nz, int ncmt, char *cmt[], int fd); +int write_2ddata(int nx, int ny, int *ntot, IMTYPE *data, int fd); +int write_3dend(int *ntot, int fd); + +#endif /* _INCLUDED_burntool_ */ Index: trunk/ppImage/src/ppImage.h =================================================================== --- trunk/ppImage/src/ppImage.h (revision 35612) +++ trunk/ppImage/src/ppImage.h (revision 35685) @@ -26,4 +26,5 @@ bool doMaskBuild; // Build internal mask bool doVarianceBuild; // Build internal variance map + bool doApplyBurntool; // apply burntool correction bool doMaskBurntool; // mask potential burntool trails bool doMaskSat; // mask saturated pixels @@ -182,5 +183,8 @@ bool ppImageCheckNoise(pmConfig *config, ppImageOptions *options, pmFPAview *view); -bool ppImageBurntoolMask(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *mask); +bool ppImageBurntoolMask(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *readout); +bool ppImageBurntoolMaskFromTable(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *readout); + +bool ppImageBurntoolApply(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *readout); bool ppImageDetrendPatternApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options); Index: trunk/ppImage/src/ppImageArguments.c =================================================================== --- trunk/ppImage/src/ppImageArguments.c (revision 35612) +++ trunk/ppImage/src/ppImageArguments.c (revision 35685) @@ -117,4 +117,10 @@ pmConfigFileSetsMD (config->arguments, &argc, argv, "LINEARITY", "-linearity", "-linearlist"); + if ((argnum = psArgumentGet(argc, argv, "-burntool"))) { + psArgumentRemove(argnum, &argc, argv); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "BURNTOOL.TABLE", PS_META_REPLACE, "", argv[argnum]); + psArgumentRemove(argnum, &argc, argv); + } + // chip selection is used to limit chips to be processed if ((argnum = psArgumentGet (argc, argv, "-chip"))) { Index: trunk/ppImage/src/ppImageBurntoolApply.c =================================================================== --- trunk/ppImage/src/ppImageBurntoolApply.c (revision 35685) +++ trunk/ppImage/src/ppImageBurntoolApply.c (revision 35685) @@ -0,0 +1,557 @@ +/* + * ppImageBurntoolApply.c + * Includes code ported from gpcsw/gpcsrc to read burntool table files and apply the fits + */ + +#include "ppImage.h" + +// Now include the files copied from gpcsw/gpcsrc + +// these get redefined +#undef MAX +#undef MIN + +#include "burntool.h" +#define EXTERN /* define EXTERN to declare variables in params.h */ +#include "burnparams.h" + + +static void * getBurntoolDataForCell(pmConfig *config, ppImageOptions *options, int burntool_cell); + +bool ppImageBurntoolApply(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *readout) +{ + // convert from pmFPA's cell number to burntool's notion + int burntool_cell = (view->cell % 8) * 8 + (view->cell - (view->cell % 8)) / 8; + + psString cellName = psMetadataLookupStr(NULL, readout->parent->concepts, "CELL.NAME"); + (void) cellName; + + CELL *burntoolData = getBurntoolDataForCell(config, options, burntool_cell); + if (!burntoolData) { + psError(PS_ERR_UNKNOWN, false, "Unable to find burntool data for cell"); + return false; + } + + psImage *image = readout->image; + +#ifdef DUMP_CELL_IMAGES + psString dumpPath = NULL; + psStringAppend(&dumpPath, "cellimages/%s.before.fits", cellName); + psphotSaveImage(NULL, image, dumpPath); + psFree(dumpPath); +#endif + + // The following code is adapted from burntool's functions burn_apply() which is in burn_fix.c and + // sub_fit() which was located in trailfit.c + + int ny = image->numRows; + + psF32 satVal = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.SATURATION"); + + for (int p = 0; p < burntoolData->npersist; p++) { + OBJBOX *box = &(burntoolData->persist[p]); + + if(box->fiterr) continue; + if(box->func != BURN_PWR && box->func != BURN_EXP) { + continue; + } + + int i, j, k; + int y0, y1, ys, dy; + + y0 = box->y0; + y1 = (box->up) ? ny-1 : 0; + dy = (y0up) { + ys = box->sy; + } else { + ys = box->ey; + } + + /* Calculate all needed y values */ + double ybuf[MAXSIZE]; + for(j=ys; dy*j<=dy*y1; j+=dy) { + if(box->func == BURN_PWR) { + if(dy*(j-y0) <= 0) { + ybuf[j] = exp(box->slope*log(Y_SCALE*1.0)); + } else { + ybuf[j] = exp(box->slope*log(Y_SCALE*dy*(j-y0))); + } + } else if(box->func == BURN_EXP) { + // ybuf[j] = exp(box->slope*Y_SCALE*MAX(0,dy*(j-y0))); + ybuf[j] = exp(box->slope*Y_SCALE*dy*(j-y0)); + } + } + + + /* Subtract the fits */ + // printf("%3d %3d %3d %1d %1d %8.5f\n", + // box->cx, box->cy, box->y0, box->func, box->up, box->slope); + + for(k=0; knfit; k++) { + i = box->xfit[k]; + for(j=box->yfit[k]; dy*j<=dy*y1; j+=dy) { + psF32 inputPixel = image->data.F32[j][i]; + if (false && inputPixel >= satVal) { + // don't modify saturated pixels + continue; + } + +#ifndef MATCH_BURNTOOL + psF32 delta = box->zero[k] * ybuf[j]; + image->data.F32[j][i] -= delta; +#else + // do math in 16 bit signed int space in order to match how the burntool program + // behaves with regard to underflows. burntool allows underflows so that the burntool + // application is reversible. + // Verified that the image matches the burntool implementation with this code + // which shows that we are reading and processing the table correctly. + int delta = box->zero[k] * ybuf[j] + 0.5; + psS16 pixel; + // adjust for BZERO + if (inputPixel >= 65535) { + pixel = 32767; + } else if (inputPixel <= 0) { + pixel = -32768; + } else { + pixel = inputPixel - 32768; + } + pixel -= delta; + // replace the pixel with correction applied + image->data.F32[j][i] = (psF32) pixel + 32768; +#endif + } + } + } +#ifdef DUMP_CELL_IMAGES + dumpPath = NULL; + psStringAppend(&dumpPath, "cellimages/%s.after.fits", cellName); + psphotSaveImage(NULL, image, dumpPath); + psFree(dumpPath); +#endif + + return true; +} + + +bool burntoolFileRead = false; +static CELL OTA[MAXCELL]; // cell structure for entire OTA +static void * getBurntoolDataForCell(pmConfig *config, ppImageOptions *options, int burntool_cell) +{ + if (!burntoolFileRead) { + psString burntoolTablePath = psMetadataLookupStr(NULL, config->arguments, "BURNTOOL.TABLE"); + if (!burntoolTablePath) { + psError(PS_ERR_UNKNOWN, false, "failed to lookup BURNTOOL.TABLE in arguments"); + return NULL; + } + psString burntoolFileName = pmConfigConvertFilename(burntoolTablePath, config, false, false); + if (!burntoolFileName) { + psError(PS_ERR_UNKNOWN, false, "failed to convert filename for %s", burntoolTablePath); + return NULL; + } + if (persist_read(OTA, burntoolFileName, 1)) { + psError(PS_ERR_UNKNOWN, "true", "failed to read burntool file"); + return NULL; + } + psFree(burntoolFileName); + burntoolFileRead = true; + } + + return (void *) &OTA[burntool_cell]; +} +#define PPIMAGE_BURNTOOL_DEBUG 0 + +// adapted from ppImageBurntoolMask +bool ppImageBurntoolMaskFromTable(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *mask) +{ + bool status = true; + int burntool_cell; + + /* Redirects and Memory juggling. */ + view->readout = 0; + psImage *image = mask->mask; + + /* Set the maskValue from the recipes. */ + psImageMaskType maskValue = options->burntoolMask; + + burntool_cell = (view->cell % 8) * 8 + (view->cell - (view->cell % 8)) / 8; + CELL *burntoolData = getBurntoolDataForCell(config, options, burntool_cell); + if (!burntoolData) { + psError(PS_ERR_UNKNOWN, false, "Unable to find burntool data for cell"); + return false; + } + +#if PPIMAGE_BURNTOOL_DEBUG + psLogMsg("ppImageBurntoolMask", 4, "Status: %d %d\n",burntoolData->npersist,maskValue); +#endif + + psLogMsg("ppImageBurntoolMask", 4, "Cell mapping: %d %d %d\n",view->cell,burntool_cell,-1); + for (int row = 0; row < burntoolData->npersist; row++) { + OBJBOX *box = &(burntoolData->persist[row]); + if (((options->burntoolTrails & 0x01)&& (box->func == 4))|| + (((options->burntoolTrails & 0x02)&& (box->up == 1))|| + ((options->burntoolTrails & 0x04)&& (box->up == 0)))) { + /* If the fit fails, burntool reports zero here. This + signifies that it expected to see a trail (else why + fit) but did not find it when it attempted to + correct. */ +#if PPIMAGE_BURNTOOL_DEBUG + psLogMsg ("ppImageBurntoolMask", 4, "Masking! %d (%d %d %d) %d %d", + burntool_cell, + ((options->burntoolTrails & 0x0001)&& (box->nfit == 0)), + ((options->burntoolTrails & 0x02)&& (box->up == 1)), + ((options->burntoolTrails & 0x04)&& (box->up == 0)), + options->burntoolTrails, + maskValue + ); +#endif + for (int i = box->sxfit; i<= box->exfit; i++) { + + if (box->up == 0) { + for (int j = 0; j <= box->y1p; j++) { +#if PPIMAGE_BURNTOOL_DEBUG + psLogMsg("ppImageBurntoolMask", 4, "Noisy!: %d %d %d %d\n", + i,j,image->data.PS_TYPE_IMAGE_MASK_DATA[j][i],maskValue); +#endif + image->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= maskValue; + } + } else { + for (int j = box->y1m; j < image->numRows ; j++) { +#if PPIMAGE_BURNTOOL_DEBUG + psLogMsg("ppImageBurntoolMask", 4, "Noisy!: %d %d %d %d\n", + i,j,image->data.PS_TYPE_IMAGE_MASK_DATA[j][i],maskValue); +#endif + image->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= maskValue; + } + } + } + } + } + + return(status); +} + +/* From persistio.c - read and write persistence info */ +/****************************************************************/ +/* persist_read(): Read all the persistence trails from a file */ +STATIC int persist_read(CELL *cell, const char *infile, int apply) +{ + int i, k, nbox=0; + char line[1024]; + FILE *fp; + +/* Initialize the counts */ + 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)); + if(boxbuf[nbox].zero == NULL || + boxbuf[nbox].xfit == NULL || + boxbuf[nbox].yfit == NULL) { + fprintf(stderr, "\rerror: failed to alloc boxbuf\n"); + exit(-673); + } + + for(i=0; i= MAXCELL) { + fprintf(stderr, "\rerror: illegal cell %d at %d from '%s'\n", + k, nbox, line); + fclose(fp); + return(-1); + } + cell[k].npersist += 1; + + if(++nbox >= MAXBURN) { + fprintf(stderr, "\rerror: persist_read overflowed MAXBURN\n"); + fclose(fp); + return(-1); + } + } + +/* Allocate some space for them */ + for(k=0; k 0) { + if( (cell[k].persist = (OBJBOX *)calloc(i, sizeof(OBJBOX))) == NULL) { + fprintf(stderr, "\rerror: failed to alloc cell persist buffer\n"); + exit(-674); + } + cell[k].npersist = 0; + } + } +/* Copy the results to the cells */ + for(i=0; inpersist; + box = cell->persist; + +/* Assess the clustering of all the streaks */ + for(i=0; i= MAXSIZE-1) continue; + for(i=box[k].sxfit; i<=box[k].exfit; i++) xusage[i] += 1; + } + +/* Identify clusters */ + for(xs=0; xs xe || box[k].exfit < xs) continue; + boxid[n] = k; + yctr[n] = (box[k].y0m+box[k].y1m+box[k].y0p+box[k].y0p+2) / 4; + n++; + } + +/* Case 1: different y start => sever */ + for(kp=0; kp 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 */ + if(zk > zj) { + if(box[j].sxfit >= box[k].sxfit) + box[j].sxfit = MAX(box[j].sxfit, box[k].exfit+1); + if(box[j].exfit <= box[k].exfit) + box[j].exfit = MIN(box[j].exfit, box[k].sxfit-1); + } else { + if(box[k].sxfit >= box[j].sxfit) + box[k].sxfit = MAX(box[k].sxfit, box[j].exfit+1); + if(box[k].exfit <= box[j].exfit) + box[k].exfit = MIN(box[k].exfit, box[j].sxfit-1); + } + } + } + } + +/* Case 2: pretty much the same y start => union */ + for(kp=0; kp exfit) */ + for(k=0; k box[k].exfit) { + for(j=k; jnpersist = nbox; + + + return(0); +} + + +/****************************************************************/ +/* persist_write(): Write all the persistence data for the next image */ +STATIC int persist_write(CELL *cell, const char *outfile) +{ + int i, k, j, err=0; + FILE *fp; + + if(outfile == NULL) return(0); + if( (fp=fopen(outfile, "w")) == NULL) { + fprintf(stderr, "\rerror: cannot open '%s' for writing\n", outfile); + return(-1); + } + +/* Dump out all the burns for the next image... */ +/* FIXME: what to do about generic cell info? */ + fprintf(fp, "# Cell: %d sky= %d rms= %d bias= %d\n", + cell[0].cell, cell[0].sky, cell[0].rms, cell[0].bias); + fprintf(fp, "#Cell time cx cy max y0 sx sy ex ey y0m y0p y1m y1p x0m x0p x1m x1p F up slope nfit sxf exf\n"); + + for(j=0; j EXPIRE_TRAIL_TIME) + continue; + + if(PERSIST_RETAIN) { +/* Keep fits which have a dubious slope */ + if(cell[j].persist[k].fiterr != FIT_SLOPE_ERROR) { + if(cell[j].persist[k].fiterr) continue; + if(cell[j].persist[k].nfit <= 0) continue; + } + } else { + if(cell[j].persist[k].fiterr) continue; + if(cell[j].persist[k].nfit <= 0) continue; + } + fprintf(fp, "%3d %7d %3d %3d %5d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %1d %1d %9.6f %3d %3d %3d %d\n", + j, cell[j].persist[k].time, + cell[j].persist[k].cx, cell[j].persist[k].cy, + cell[j].persist[k].max, cell[j].persist[k].y0, + cell[j].persist[k].sx, cell[j].persist[k].sy, + cell[j].persist[k].ex, cell[j].persist[k].ey, + cell[j].persist[k].y0m, cell[j].persist[k].y0p, + cell[j].persist[k].y1m, cell[j].persist[k].y1p, + cell[j].persist[k].x0m, cell[j].persist[k].x0p, + cell[j].persist[k].x1m, cell[j].persist[k].x1p, + cell[j].persist[k].func, cell[j].persist[k].up, + cell[j].persist[k].slope, cell[j].persist[k].nfit, + cell[j].persist[k].sxfit, cell[j].persist[k].exfit, + cell[j].persist[k].fiterr); + for(i=0; idoMaskBuild) { @@ -74,6 +74,19 @@ pmMaskBadPixels(input, mask, options->maskValue); } + + if (options->doApplyBurntool) + { + // extern bool ppImageBurntoolApply(pmConfig *, ppImageOptions *, pmFPAview *, pmReadout *); + ppImageBurntoolApply(config, options, view, input); + } + if (options->doMaskBurntool) { - ppImageBurntoolMask(config,options,view,input); + if (options->doApplyBurntool) { + // build burntool mask from data input burntool table + ppImageBurntoolMaskFromTable(config,options,view,input); + } else { + // build burntool mask from data in the input image's fits extension + ppImageBurntoolMask(config,options,view,input); + } } Index: trunk/ppImage/src/ppImageOptions.c =================================================================== --- trunk/ppImage/src/ppImageOptions.c (revision 35612) +++ trunk/ppImage/src/ppImageOptions.c (revision 35685) @@ -23,4 +23,5 @@ options->doMaskLow = false; // mask low pixels options->doMaskBurntool = false; // mask potential burntool trails + options->doApplyBurntool = false; // apply burntool correction options->doVarianceBuild = false; // Build internal variance options->doMask = false; // Mask bad pixels @@ -231,4 +232,5 @@ options->doMaskLow = psMetadataLookupBool(NULL, recipe, "MASK.LOW"); options->doMaskBurntool = psMetadataLookupBool(NULL, recipe, "MASK.BURNTOOL"); + options->doApplyBurntool = psMetadataLookupBool(NULL, recipe, "APPLY.BURNTOOL"); options->doVarianceBuild = psMetadataLookupBool(NULL, recipe, "VARIANCE.BUILD"); options->doAuxMask = psMetadataLookupBool(NULL, recipe, "MASK.AUXMASK"); Index: trunk/ppImage/src/ppImageParseCamera.c =================================================================== --- trunk/ppImage/src/ppImageParseCamera.c (revision 35612) +++ trunk/ppImage/src/ppImageParseCamera.c (revision 35685) @@ -184,4 +184,12 @@ } outCTE->save = true; + } + if (options->doApplyBurntool) { + // If no burntool file was supplied do not fail. (camera_exp.pl does not supply it + // for example. + if (!psMetadataLookupStr(NULL, config->arguments, "BURNTOOL.TABLE")) { + psWarning("BURNTOOL.TABLE not supplied setting doApplyBurntool to false"); + options->doApplyBurntool = false; + } }