Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 42723) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 42724) @@ -17,18 +17,16 @@ /* burn_fix(): Find and fix all the burned patches in a cell */ STATIC int burn_fix(int nx, int ny, int NX, int NY, IMTYPE *buf, - CELL *cell, int cellnum) + CELL *cell, int cellnum, char *camera) { int k, err; /* Find all the burned patches and bright stars */ - err = star_detect(nx, ny, NX, NY, imbuf, mbuf, msbuf, cell, cellnum); - - err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell); + err = star_detect(nx, ny, NX, NY, imbuf, mbuf, msbuf, cell, cellnum, camera); + + err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell, camera); // 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); -// err = grow_mask(nx, ny, NX, mbuf, XMASK_GROW, YMASK_GROW, cell->nburn, cell->burn); err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW, MASK_SAT_HALO, cell->nburn, cell->burn); @@ -44,5 +42,5 @@ // 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); + cell->sky+cell->bias, cell->rms, BURN_PWR, camera); /* Subtract the fit from the image */ // fprintf(stderr, "Subtracting trail %d\n", k); @@ -54,4 +52,5 @@ for(k=0; krms, mask, cell->burn+k); + err = burn_test(nx, ny, NX, data, cell->rms, mask, cell->burn+k, camera); } @@ -240,8 +239,8 @@ */ STATIC int burn_test(int nx, int ny, int NX, DTYPE *data, int rms, - MTYPE *mask, OBJBOX *box) + MTYPE *mask, OBJBOX *box, char *camera) { int i, j, nburn, nref, nmed; - int y0, y1, ymid, xburn, xref, dx=0, dy=0; + int y0, y1, x0, x1, xwidth, ymid, xburn, xref, dx=0, dy=0; /*introduce a maximum to the burn contrast check, to avoid a very large number in crowded regions*/ @@ -275,4 +274,6 @@ y0 = ymid + MAX(box->ey-ymid, FIT_EDGE); y1 = MIN(y0+100, ny-1); + x0 = box->sx; + x1 = box->ex; if(2*box->ex - box->sx + 1 < nx-1) { /* Work right? */ @@ -294,8 +295,28 @@ dx = 0; dy = MAX(ymid-box->sy, box->ey-ymid); + + if(!strcmp(camera,"gpc2")) { + /*start right at the edge of the star box might not be the best*/ + /*add 20% of box size on top, with a maximum of ten pixels*/ + dy = dy + MIN((box->ey-box->sy)/5,10); + + /*Use the inner third of the star width for this test. Should be brightest there*/ + xwidth = MAX((box->ex-box->sx)/3,10); + x0 = box->sx + (box->ex-box->sx)/2 - xwidth/2; + x1 = box->sx + (box->ex-box->sx)/2 + xwidth/2; + } else { + //regular gpc1 setup + x0 = box->sx; + x1 = box->ex; + } + if(dy < FIT_EDGE) dy = FIT_EDGE; y0 = ymid + dy; y1 = y0 + MIN(ny-1-y0, 2*ymid-y0-1); if(y1-y0 > test_box) y1 = y0 + test_box; + + } + if(VERBOSE & VERB_BOXGROW) { + printf("test input dx,dy,x0,x1, y0, y1, test_box, ymid: %8d %8d %8d %8d %8d %8d %8d %8d\n", dx,dy,x0,x1, y0, y1, test_box, ymid); } @@ -305,5 +326,5 @@ xburn = xref = 0; nburn = nref = 0; - for(i=box->sx; i<=box->ex; i++) { + for(i=x0; i<=x1; i++) { if(mask[i+j*NX] == MASK_NONE) { xburn += data[i+j*NX]; @@ -326,5 +347,5 @@ if(nref > 0 && nburn > 0) median_buf[nmed++] = xburn - xref; if(VERBOSE & VERB_BOXGROW) { - printf("contrast: %8d %8d %8d %8d %8d\n", j, xburn, (2*ymid-j), xref, xburn - xref); + printf("contrast: %8d %8d %8d %8d %8d %8d %8d\n", j, xburn,nburn, (2*ymid-j), xref,nref, xburn - xref); } @@ -332,4 +353,8 @@ box->diff = int_median(nmed, median_buf); box->burned = (box->diff > burn_contrast); + if(VERBOSE & VERB_BOXGROW) { + printf("test output: %8d %8d %8d\n", box->burned, burn_contrast, box->diff); + } + return(0); } @@ -344,5 +369,5 @@ /* 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) + MTYPE *mask, OBJBOX *box, char *camera) { int i, j, n, nrow, nmed; Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 42723) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 42724) @@ -29,5 +29,5 @@ HeaderUnit ehu; const char* ifilename = "-"; - char extname[FH_MAX_STRLEN+1], otaposn[FH_MAX_STRLEN+1]; + char extname[FH_MAX_STRLEN+1], otaposn[FH_MAX_STRLEN+1], *camera="gpc1"; int cellmask[MAXCELL]; int otanum; @@ -176,4 +176,8 @@ } else if(strncmp(argv[i], "infits=", 7) == 0) { /* infits=fname */ persistfitsfile = argv[i] + 7; + +/* Input camera keyword (optional) */ + } else if(strncmp(argv[i], "camera=", 7) == 0) { /* camera=fname */ + camera = argv[i] + 7; /* Keep persistence streaks which had a bad slope? */ @@ -577,5 +581,5 @@ /* Use the table-driven fits instead of calculating new ones */ burn_apply(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, - buf, OTA+cell); + buf, OTA+cell,camera); /* Tell us about it? */ if(VERBOSE & VERB_NORM) burn_blab(OTA+cell); @@ -584,10 +588,10 @@ /* Restore the old burns */ burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, - buf, OTA+cell); + buf, OTA+cell,camera); } else { /* Fix up the burns */ burn_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, naxis2, buf, - OTA+cell, cell); + OTA+cell, cell,camera); /* Collect up a good star list */ @@ -602,5 +606,5 @@ persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf, - OTA+cell); + OTA+cell,camera); /* Tell us about it? */ if(VERBOSE & VERB_NORM) persist_blab(OTA+cell); Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 42723) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 42724) @@ -172,12 +172,12 @@ 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); + CELL *cell, int cellnum,char *camera); 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); + MTYPE *mask, OBJBOX *box,char *camera); +STATIC int burn_restore(int nx, int ny, int NX, IMTYPE *buf, CELL *cell,char *camera); +STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell,char *camera); 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_fix(int nx, int ny, int stride, IMTYPE *buf, CELL *cell,char *camera); STATIC int persist_merge(CELL *cell); @@ -187,7 +187,7 @@ STATIC int star_detect(int nx, int ny, int NX, int NY, DTYPE *data, - MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum); + MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum,char *camera); STATIC int burn_check(int nx, int ny, int stride, int NY, DTYPE *buf, - MTYPE *mask, CELL *cell); + MTYPE *mask, CELL *cell,char *camera); STATIC int cell_stats(int nx, int ny, int NX, int NY, DTYPE *data, CELL *cell); @@ -197,5 +197,5 @@ STATIC int fit_trail(int nx, int ny, int NX, DTYPE *data, MTYPE *mask, - OBJBOX *box, int up, int sky, int rms, int fitfunc); + OBJBOX *box, int up, int sky, int rms, int fitfunc,char *camera); STATIC int sub_fit(int nx, int ny, int NX, IMTYPE *buf, OBJBOX *box, int sign); Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c (revision 42723) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c (revision 42724) @@ -18,5 +18,5 @@ /* burn_fix must have been run first! */ STATIC int persist_fix(int nx, int ny, int NX, IMTYPE *buf, - CELL *cell) + CELL *cell, char *camera) { int k, err=0; @@ -43,5 +43,5 @@ /* Fit the trail */ err = fit_trail(nx, ny, NX, imbuf, mbuf, cell->persist+k, 0, - cell->sky+cell->bias, cell->rms, BURN_EXP); + cell->sky+cell->bias, cell->rms, BURN_EXP, camera); /* Subtract out the fit */ Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c (revision 42723) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c (revision 42724) @@ -17,5 +17,5 @@ /* star_detect(): Find all the stars and burned patches in a cell */ STATIC int star_detect(int nx, int ny, int NX, int NY, DTYPE *data, - MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum) + MTYPE *mask, MTYPE *veto, CELL *cell, int cellnum, char *camera) { int i, j, k, l; @@ -124,5 +124,5 @@ /* But take a hard look at stars with really bright centers... */ if(!boxbuf[nbox].burned && boxbuf[nbox].max > MAX_THRESH) { - burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox); + burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox, camera); } boxbuf[nbox].time = cell->time; Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c (revision 42723) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c (revision 42724) @@ -14,5 +14,5 @@ #include "burnparams.h" -STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE]; +STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE],xbuf[MAXSIZE]; /****************************************************************/ @@ -92,9 +92,9 @@ STATIC int fit_trail(int nx, int ny, int NX, DTYPE *data, MTYPE *mask, OBJBOX *box, - int up, int bckgnd, int rms, int fitfunc) + int up, int bckgnd, int rms, int fitfunc,char *camera) { - int i, j, k, err, nfit,firstfitpix; - int xs, xe, y0, y1, y2, dy, yfit; - double slope, zero, zsum, wsum, trial,zfun; + int i, j, k, err, nfit,firstfitpix,secondfitpix, xwidth; + int xs, xe, y0, y1, y2, dy, yfit,x0,x1; + double slope, zero, zsum, wsum, trial,zfun,xsum,xmean,xstd; char fooname[60]; FILE *fp=stdout; @@ -104,5 +104,5 @@ return(-1); } - + /*xs = box->sxfit; */ /*xe = box->exfit; */ @@ -110,4 +110,5 @@ xs = box->sx; xe = box->ex; + if(up) { @@ -119,4 +120,5 @@ yfit = box->sy; /* Start of validity, TBD later by col */ firstfitpix = MAX((box->ey - box->sy)*2,20) ; /* Number of pixels to include in the initial fit*/ + secondfitpix = 200; /* Number of pixels to include in the second fit*/ } else { //y0 = (box->y0m + box->y1m + box->y0p + box->y1p + 2) / 4; @@ -127,10 +129,24 @@ yfit = box->ey; /* Start of validity, TBD later by col */ firstfitpix = 100; /* Number of pixels to include in the initial fit*/ - + secondfitpix = 200; /* Number of pixels to include in the second fit*/ + } + + if(!strcmp(camera,"gpc2")) { + /*Use the inner third of the star width for this test. Should be brightest there*/ + xwidth = MAX((box->ex-box->sx)/3,10); + x0 = box->sx + (box->ex-box->sx)/2 - xwidth/2; + x1 = box->sx + (box->ex-box->sx)/2 + xwidth/2; + /*when fitting a burn, do not want to include part of the star*/ + /*add 20% of box size on top, with a maximum of ten pixels*/ + if(up) y1 = box->ey + MIN((box->ey-box->sy)/5,10); + } else { + //regular gpc1 setup + x0 = box->sx; + x1 = box->ex; } if(VERBOSE & VERB_FIT) { - printf("Fit requested xs,xe,y0,y1,y2,dy,up,f= %d %d %d %d %d %d %d %d\n", - xs, xe, y0, y1, y2, dy, up, fitfunc); + printf("Fit requested xs,xe,y0,y1,y2,dy,up,f,x0,x1= %d %d %d %d %d %d %d %d %d %d\n", + xs, xe, y0, y1, y2, dy, up, fitfunc,x0,x1); } @@ -184,5 +200,5 @@ ybuf[j] = Y_SCALE * dy * (j-y0); zbuf[j]=0.0; - for(i=xs, k=0; i<=xe; i++) { + for(i=x0, k=0; i<=x1; i++) { if(data[i+j*NX] != NODATA && (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO) ) { zbuf[j] += data[i+j*NX]; @@ -221,9 +237,9 @@ /* Trim away any really big positive deviations */ for(j=y1; j<=y2; j++) { -// if(zbuf[j]-(ybuf[j]*slope+zero) > MAXDEV && -// zbuf[j] > log(MAXRMS*rms)) wbuf[j] = 0; if(abs(zbuf[j]-(ybuf[j]*slope+zero)) > MAXDEV) wbuf[j] = 0; - if(VERBOSE & VERB_FITPROF) { - fprintf(fp, "%9.4f %9.4f %9.4f %3d\n", ybuf[j], zbuf[j], wbuf[j], j); + if(j < (y2-secondfitpix) && !up) wbuf[j] = 0; + if(j > (y1+secondfitpix) && up) wbuf[j] = 0; + if(VERBOSE & VERB_FIT) { + printf("trim j : %d %9.4f %9.4f %9.4f %d %9.4f %9.4f\n", j,ybuf[j],zbuf[j],(ybuf[j]*slope+zero),abs(zbuf[j]-(ybuf[j]*slope+zero)), wbuf[j], wbuf0[j]); } } @@ -237,5 +253,5 @@ if(VERBOSE & VERB_FITPROF) { - printf("jt %d %.3f %.3f\n", xs, zero, slope); + printf("2nd fit nfit, zero, slope: %d %.3f %.3f\n", nfit, zero, slope); fclose(fp); } @@ -249,4 +265,5 @@ linearrms(nfit, ybuf+y1, zbuf+y1, slope, zero, &trial); /* 100203 JT: bad idea: appears to be a bug in read/writing ABS(nfit) */ + // if(trial > 2*rms) box->nfit = -box->nfit; if(trial > 2*rms) box->func = BURN_POSSLOPE; @@ -256,50 +273,98 @@ } -/* Reconstruct the fit for comparison with the data col by col */ - for(j=y1; j<=y2; j++) zbuf[j] = exp(ybuf[j]*slope); - /* And tag on the extras to determine the end of the fit */ - for(j=yfit; dy*j NEGLIGIBLE_TRAIL*rms) break; + if(zfun < NEGLIGIBLE_TRAIL*rms) break; + } } //For burns, just let it go to the top - if(j < y2 && !up) box->eyfit = j; - -/* Determine good x-ranges and Save the fit information */ - box->slope = slope; - box->nfit = 0; -/* Save individual scaled versions for each column */ - for(i=xs; i<=xe; i++) { - zsum = wsum = 0.0; - for(j=y1; j<=y2; j++) { + //if(j < y2 && !up) box->eyfit = j; + box->eyfit = j; + + + if(!strcmp(camera,"gpc2")) { + + /*Conpute the mean and standard deviation of the trail along x*/ + for(i=xs; i<=xe; i++) { + xbuf[i] = 0.0; + for(j=y1, k=0; j<=y2; j++) { + if(wbuf[j] > 0 && data[i+j*NX] != NODATA && + (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO)) { + xbuf[i] += data[i+j*NX]; + k++; + } + } + xbuf[i] = xbuf[i] / MAX(1,k); + } + + xsum = wsum = 0.; + for(i=xs; i<=xe; i++) { + xsum += i*xbuf[i]; + wsum += xbuf[i]; + } + xmean = xsum / wsum; + + xsum = wsum = 0.; + for(i=xs; i<=xe; i++) { + xsum += xbuf[i]*pow(i-xmean,2); + wsum += xbuf[i]; + } + xstd = sqrt(xsum / (wsum-1.)); + if(VERBOSE & VERB_FIT) { + printf("x mean xstd : %9.4f %9.4f %9.4f %9.4f\n", xmean, xstd,xsum,wsum); + } + + /* Update fit ranges */ + if(xstd > 1 && xstd < 50) { + box->sxfit = round(xmean-(xstd*1.1)); + box->exfit = round(xmean+(xstd*1.1)); + } else { + box->sxfit = xs; + box->exfit = xs - 1; + box->fiterr = FIT_ALL_GONE; + } + } else { + /* Reconstruct the fit for comparison with the data col by col */ + for(j=y1; j<=y2; j++) zbuf[j] = exp((ybuf[j]*slope)); + /* And tag on the extras to determine the end of the fit */ + for(j=yfit; dy*jslope = slope; + box->nfit = 0; + /* Save individual scaled versions for each column */ + for(i=xs; i<=xe; i++) { + zsum = wsum = 0.0; + for(j=y1; j<=y2; j++) { if(wbuf[j] > 0 && data[i+j*NX] != NODATA && (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO)) { -// if(wbuf[j] > 0 && data[i+j*NX] != NODATA) { -// wsum += wbuf[j]; -// zsum += wbuf[j] * (data[i+j*NX] - bckgnd) / zbuf[j]; wsum += zbuf[j]; zsum += data[i+j*NX] - bckgnd; } - } - if(zsum < 0 || wsum <= 0) { + } + if(zsum < 0 || wsum <= 0) { zsum = 0.0; - } else { + } else { zsum = zsum / wsum; - } - /* FIXME: what's a really good criterion for negligible fit? */ - /* 100 pixels up fit is zsum or ~zsum/e */ -// if(zsum > NEGLIGIBLE_TRAIL*rms) { - /* 100113: but also evaluate at the end of the fit for stubby trails */ - if(zsum > NEGLIGIBLE_TRAIL*rms || + } + /* FIXME: what's a really good criterion for negligible fit? */ + /* 100 pixels up fit is zsum or ~zsum/e */ + /* 100113: but also evaluate at the end of the fit for stubby trails */ + //printf("scaled fit comparison x zsum zbuf[y2]*zsum thresh: %d %9.4f %9.4f %9.4f %f %d %f\n", i, zsum, zbuf[y2], zbuf[y2]*zsum, NEGLIGIBLE_TRAIL, rms, NEGLIGIBLE_TRAIL*rms); + + if(zsum > NEGLIGIBLE_TRAIL*rms || zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) { /* Ascertain the starting point of where the fit is good */ @@ -314,15 +379,24 @@ box->zero[box->nfit] = zsum; box->xfit[box->nfit] = i; + //printf("y comparison : %d %d %9.4f %9.4f %d %d\n", i, box->yfit[box->nfit], box->zero[box->nfit],zsum, box->xfit[box->nfit],box->nfit); box->nfit += 1; - } - } -/* Update fit ranges */ - 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; + } + } + /* Update fit ranges */ + 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(VERBOSE & VERB_FIT) { + printf("Final fit params xsfit,xefit,y1,y2, yfit,slope,nfit = %d %d %d %d %d %.3f %d\n", + box->sxfit, box->exfit, y1, y2,box->eyfit, box->slope, box->nfit); }