Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/Makefile =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/Makefile (revision 24302) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/Makefile (revision 24391) @@ -32,5 +32,5 @@ pscoords/pscoords.h -$(OBJ)/psfstats.o: psfstats.c burntool.h +$(OBJ)/psfstats.o: psfstats.c burntool.h psf/psf.h $(OBJ)/sort.o: sort.c Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 24302) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 24391) @@ -105,4 +105,19 @@ if(!cell->persist[k].fiterr) sub_fit(nx, ny, NX, buf, &(cell->persist[k]), -1); + } + return(0); +} + + +/****************************************************************/ +/* burn_apply(): Subtract the trail fits from the image */ +STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell) +{ + int k; + +/* Restore all the burns */ + for(k=0; knpersist; k++) { + if(!cell->persist[k].fiterr) + sub_fit(nx, ny, NX, buf, &(cell->persist[k]), +1); } return(0); Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 24302) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 24391) @@ -33,5 +33,5 @@ int otanum; int nextend, cellxy, cell, cellcode; - int ext, update, restore, psfsize, psfavg; + int ext, update, restore, apply, tableonly, psfsize, psfavg; IMTYPE *buf; const char *burnfile=NULL, *persistfile=NULL, *persistfitsfile=NULL; @@ -102,6 +102,8 @@ /* Parse the args */ cellxy = -1; - update = 1; - restore = 0; + update = 1; /* Calc fits, apply fits, write img and table */ + restore = 0; /* Restore previous fit only, write img */ + apply = 0; /* Apply previous fit only, write img */ + tableonly = 0; /* Calc fits (apply fits), write table only */ psfsize = 32; psfavg = 0; @@ -137,7 +139,15 @@ update = argv[i][7] == 'y' || argv[i][7] == '1' || argv[i][7] == 't'; -/* Modify the input MEF by subtracting fits? */ +/* Modify the input MEF by adding back fits? */ } else if(strncmp(argv[i], "restore=", 8) == 0) { /* restore={t|f} */ restore = argv[i][8] == 'y' || argv[i][8] == '1' || argv[i][8] == 't'; + +/* Modify the input MEF by subtracting previously calculated fits? */ + } else if(strncmp(argv[i], "apply=", 6) == 0) { /* apply={t|f} */ + apply = argv[i][6] == 'y' || argv[i][6] == '1' || argv[i][6] == 't'; + +/* Calculate and write tables only? */ + } else if(strncmp(argv[i], "tableonly=", 10) == 0) {/* tableonly={t|f} */ + tableonly = argv[i][10] == 'y' || argv[i][10] == '1' || argv[i][10] == 't'; /* Output file for burn streaks */ @@ -150,9 +160,9 @@ /* Input text file for previous burn persistence streaks */ - } else if(strncmp(argv[i], "infits=", 8) == 0) { /* in=fname */ + } else if(strncmp(argv[i], "trailin=", 8) == 0) {/* trailin=fname */ persistfile = argv[i] + 8; /* Same thing, but information is stored in tables in a FITS file. */ - } else if(strncmp(argv[i], "trailinfits=", 8) == 0) { /* infits=fname */ + } else if(strncmp(argv[i], "trailinfits=", 12) == 0) { /* trailinfits=fname */ persistfitsfile = argv[i] + 12; @@ -162,6 +172,6 @@ /* Same thing, but information is stored in tables in a FITS file. */ - } else if(strncmp(argv[i], "trailinfits=", 8) == 0) { /* trailin=fname */ - persistfitsfile = argv[i] + 12; + } else if(strncmp(argv[i], "infits=", 7) == 0) { /* infits=fname */ + persistfitsfile = argv[i] + 7; /* Output file for PSF gallery */ @@ -322,6 +332,6 @@ /* If there is no other persistence info supplied, try getting * it from the input FITS file. */ - if(restore && persistfile == NULL) { - if(persistfitsfile == NULL) persistfitsfile = ifilename; + if((restore || apply) && persistfile == NULL) { + if(persistfitsfile == NULL) persistfitsfile = ifilename; } @@ -334,9 +344,9 @@ /* Read the persistence data for this OTA */ - if(persistfile != NULL) { + if(persistfile != NULL) { /* Text data file */ if(persist_read(OTA, persistfile)) exit(EXIT_FAILURE); - } - else if(persistfitsfile != NULL) { - if(persist_fits_read(OTA, persistfitsfile) != FH_SUCCESS) exit(EXIT_FAILURE); + } else if(persistfitsfile != NULL) { /* FITS table */ + if(persist_fits_read(OTA, persistfitsfile) != FH_SUCCESS) + exit(EXIT_FAILURE); } @@ -369,5 +379,5 @@ "warning: Restoring old burns, but header indicates no burns previously corrected.\n"); } - else if (update && (burn_applied == FH_TRUE)) { + else if ((update||apply) && (burn_applied == FH_TRUE)) { fprintf(stderr, "warning: Applying burn correction, but header indicates burns previously corrected.\n"); @@ -491,5 +501,5 @@ if(VERBOSE & VERB_NORM) { - printf("nx=%d ny=%d prex=%d postx=%d posty=%d BZERO=%.1f\n", + printf("nx=%d ny=%d prex=%d postx=%d posty=%d BZERO=%d\n", naxis1, naxis2, prescan1, ovrscan1, ovrscan2, BZERO); } @@ -542,7 +552,17 @@ } - - if(!restore) { - + if(apply) { +/* Use the table-driven fits instead of calculating new ones */ + burn_apply(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, + buf, OTA+cell); +/* Tell us about it? */ + if(VERBOSE & VERB_NORM) burn_blab(OTA+cell); + + } else if(restore) { +/* Restore the old burns */ + burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, + buf, OTA+cell); + + } else { /* Fix up the burns */ burn_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, naxis2, buf, @@ -559,19 +579,15 @@ if(VERBOSE & VERB_NORM) burn_blab(OTA+cell); -/* Fix up the streaks */ - persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf, - OTA+cell); - +/* Fix up the streaks (don't bother if table only) */ + if(!tableonly) { + persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf, + OTA+cell); /* Tell us about it? */ - if(VERBOSE & VERB_NORM) persist_blab(OTA+cell); - - } else { -/* Restore the old burns */ - burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, - buf, OTA+cell); + if(VERBOSE & VERB_NORM) persist_blab(OTA+cell); + } } /* Write the corrected data back to the FITS. */ - if(update) { + if(!tableonly && update) { fh_ehu(ehu, 0); /* Seek back to the start of data */ if (fh_write_padded_image(ehu, fh_file_desc(ehu), buf, @@ -590,10 +606,10 @@ /* Dump out the postage stamp file */ - if(psffile != NULL) { + if(psffile != NULL && !restore && !apply) { psf_write(psfsize, psfsize, OTA, otanum, psffile); } /* Dump out the PSF stats */ - if(psfstatfile != NULL) { + if(psfstatfile != NULL && !restore && !apply) { psf_write_stats(psfsize, psfsize, OTA, otanum, psfstatfile, psfavg); } @@ -601,11 +617,15 @@ /* Write burn info to FITS file. */ if(update) persist_fits_write(OTA, ihu); - - if(restore) { + if(restore || tableonly) { /* Indicate in the header that the burns are not applied. */ - fh_set_bool(ihu, FH_AUTO, PHU_NAME_BURN_APPLIED, - FH_FALSE, PHU_COMMENT_BURN_APPLIED); - fh_rewrite(ihu); + fh_set_bool(ihu, FH_AUTO, PHU_NAME_BURN_APPLIED, + FH_FALSE, PHU_COMMENT_BURN_APPLIED); + fh_rewrite(ihu); + } else if(apply) { + /* Indicate in the header that the burns have been applied. */ + fh_set_bool(ihu, FH_AUTO, PHU_NAME_BURN_APPLIED, + FH_TRUE, PHU_COMMENT_BURN_APPLIED); + fh_rewrite(ihu); } @@ -711,6 +731,8 @@ printf(" cell=N Work on just one cell? Cell count [0:63] mode.\n"); printf(" mask=0101... 64 digits to work on cells 0:63.\n"); - printf(" update={t|f} Modify the input MEF by subtracting fits?\n"); + printf(" update={t|f} Modify the input MEF writing table and subtracting fits?\n"); printf(" restore={t|f} Restore the input MEF by adding input fits?\n"); + printf(" apply={t|f} Modify the input MEF by subtracting previously calculated fits?\n"); + printf(" tableonly={t|f} Calculate fits but do *not* modify the input MEF images, only write tables\n"); printf(" in=fname Input file for previous burn persistence streaks\n"); printf(" infits=fname Input FITS file for previous burn persistence streaks (stored\n"); Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 24302) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 24391) @@ -173,4 +173,5 @@ 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); STATIC int persist_write(CELL *cell, const char *outfile); Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 24302) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 24391) @@ -66,7 +66,54 @@ inhibited from subtracting the fits by "update=f". + Two other modes can be invoked using "tableonly=t" and "apply=t". + The first mode is just like "update=t", except that it does not write + the modified images, only the FITS table and/or tabulated fit data. + The second mode is like "restore=t", except that it subtracts the + tabulated fits instead of adding them back in. Note that this requires + slightly different specifications for the "in=" file: + + update or tableonly: in=previous_file_for_persistence + restore or apply: in=this_file_for_burn_and_persistence + + Thus a sequence of images could be processed (leaving off usual + arguments such as mask=, etc) to immediately fit and subtract: + + burntool o4991g0001o24.fits \ + out=o4991g0001o24.burn \ + psfstat=o4991g0001o24.stat psf=o4991g0001o24.psf + + burntool o4991g0002o24.fits \ + in=o4991g0001o24.burn out=o4991g0002o24.burn \ + psfstat=o4991g0002o24.stat psf=o4991g0002o24.psf + + burntool o4991g0003o24.fits \ + in=o4991g0002o24.burn out=o4991g0003o24.burn \ + psfstat=o4991g0003o24.stat psf=o4991g0003o24.psf + + or as a two step process of generating tables: + + burntool o4991g0001o24.fits tableonly=t \ + out=o4991g0001o24.burn \ + psfstat=o4991g0001o24.stat psf=o4991g0001o24.psf + + burntool o4991g0002o24.fits tableonly=t \ + in=o4991g0001o24.burn out=o4991g0002o24.burn \ + psfstat=o4991g0002o24.stat psf=o4991g0002o24.psf + + burntool o4991g0003o24.fits tableonly=t \ + in=o4991g0002o24.burn out=o4991g0003o24.burn \ + psfstat=o4991g0003o24.stat psf=o4991g0003o24.psf + + and then later on carrying out the subtractions of the fits: + + burntool o4991g0001o24.fits apply=t + burntool o4991g0002o24.fits apply=t + burntool o4991g0003o24.fits apply=t + The identification of significant (but unsaturated) stars is an important component of burntool's mask generation, and if requested burntool can assemble a gallery of postage stamps of suitable stars. + This does not work with "restore=t" or "apply=t" of course, since + these work only on previous tables and do not examine the images. Any star whose maximum exceeds PSF_THRESH above sky but is not deemed @@ -125,7 +172,16 @@ qt=-5.201 Tangential quadrupole [pix^2] (-99.99 if fails) q+ * cos(2*phi) + qx * sin(2*phi) + q3c=-1.193 Cosine trefoil measure [pix^2] (-99.99 if fails) + q3s=-0.038 Sine trefoil measure [pix^2] (-99.99 if fails) qpm=5.154 Plus quadrupole averaged over psfavg qcm=0.560 Cross quadrupole averaged over psfavg qtm=-5.201 Tangential averaged over psfavg + q3cm=-1.230 Cosine trefoil averaged over psfavg + q3sm=-0.079 Sine trefoil averaged over psfavg + + The net, polar coordinate quadrupole can be assembled from + + Qmagnitude = sqrt(qpm^2+qcm^2) + Qtheta = 0.5 * atan2(qcm, qpm) Helpful utilities include: @@ -151,8 +207,18 @@ update={t|f} - Modify the input MEF by subtracting fits (default t)? + Calculate fits and modify the input MEF images by subtracting + the trail fits and writing a FITS table (default t). If + false the MEF is not altered, although PSF can be calculated. restore={t|f} - Restore the input MEF by adding input fits (default f)? + Restore the input MEF by adding input fits (default f) + + apply={t|f} + Use previously calculated fits to modify the input MEF image + by subtraction (default f) + + tableonly={t|f} + Calculate fits but do *not* modify the input MEF images, + only write the fit data as a FITS and/or text table (default f) in=fname Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/psfstamp.c =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/psfstamp.c (revision 24302) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/psfstamp.c (revision 24391) @@ -136,5 +136,5 @@ int i, k, l, nstar, fdout, otacx, otacy, xid, yid, ntot=0, sumax; int cellcount, ota_xid, ota_yid; - double scale, phi, fwhm[3], q[3], qt, xfp, yfp, pi=4*atan(1.0); + double scale, phi, fwhm[3], q[5], qt, xfp, yfp, pi=4*atan(1.0); IMTYPE *median_image; CELL *cell; @@ -203,6 +203,7 @@ qt = -q[1] * cos(2*phi) - q[2] * sin(2*phi); - printf("N= %d PSFmaj= %.2f min= %.2f theta= %.1f m2= %.2f q+= %.3f qx= %.3f qt= %.3f\n", - nstar, fwhm[0], fwhm[1], fwhm[2]*180/pi, q[0], q[1], q[2], qt); + printf("N= %d PSFmaj= %.2f min= %.2f theta= %.1f m2= %.2f q+= %.3f qx= %.3f qt= %.3f q3c= %.3f q3s= %.3f\n", + nstar, fwhm[0], fwhm[1], fwhm[2]*180/pi, q[0], q[1], q[2], qt, + q[3], q[4]); if(CONCAT_FITS) { @@ -258,9 +259,12 @@ double m2[MAXPSFMEDIAN], qp[MAXPSFMEDIAN], qc[MAXPSFMEDIAN]; double qt[MAXPSFMEDIAN], fwavg[MAXPSFMEDIAN]; + double q3c[MAXPSFMEDIAN], q3s[MAXPSFMEDIAN]; double qpavg[MAXPSFMEDIAN], qcavg[MAXPSFMEDIAN], qtavg[MAXPSFMEDIAN]; + double q3cavg[MAXPSFMEDIAN], q3savg[MAXPSFMEDIAN]; int nstar[MAXCELL], nfw[MAXCELL]; double fwmed[MAXCELL], m2med[MAXCELL]; double qpmed[MAXCELL], qcmed[MAXCELL], qtmed[MAXCELL]; - double qpmacro, qcmacro, qtmacro, fwmacro; + double q3cmed[MAXCELL], q3smed[MAXCELL]; + double qpmacro, qcmacro, qtmacro, fwmacro, q3cmacro, q3smacro; FILE *fp; CELL *cell; @@ -304,4 +308,6 @@ qp[nstar[k]] = q[1]; qc[nstar[k]] = q[2]; + q3c[nstar[k]] = q[3]; + q3s[nstar[k]] = q[4]; /* Get the position in the focal plane and therefore the qt statistic */ psc_cell_to_pixel(cellx, celly, 0.5*PSC_HCELL/PSC_PIXEL, @@ -325,6 +331,9 @@ qcmed[k] = double_median(nstar[k], qc); qtmed[k] = double_median(nstar[k], qt); + q3cmed[k] = double_median(nstar[k], q3c); + q3smed[k] = double_median(nstar[k], q3s); } else { m2med[k] = qpmed[k] = qcmed[k] = qtmed[k] = -99.99; + q3cmed[k] = q3smed[k] = -99.99; } /* Toss these results into the macrocell median hopper */ @@ -337,4 +346,6 @@ qcavg[nqavg] = qc[l]; qtavg[nqavg] = qt[l]; + q3cavg[nqavg] = q3c[l]; + q3savg[nqavg] = q3s[l]; nqavg++; } @@ -346,5 +357,7 @@ qcmacro = double_median(nqavg, qcavg); qtmacro = double_median(nqavg, qtavg); - if(nqavg == 0) qpmacro = qcmacro = qtmacro = -99.99; + q3cmacro = double_median(nqavg, q3cavg); + q3smacro = double_median(nqavg, q3savg); + if(nqavg == 0) qpmacro = qcmacro = qtmacro = q3cmacro = q3smacro = -99.99; fwmacro = double_median(nfwave, fwavg); @@ -357,9 +370,9 @@ cell = ota + k; - fprintf(fp, "ext=xy%1d%1d bias=%d sky=%d rmssky=%d npsf=%d fwhm=%.2f fwmed=%.2f m2=%.2f qp=%.3f qc=%.3f qt=%.3f qpm=%.3f qcm=%.3f qtm=%.3f\n", + fprintf(fp, "ext=xy%1d%1d bias=%d sky=%d rmssky=%d npsf=%d fwhm=%.2f fwmed=%.2f m2=%.2f qp=%.3f qc=%.3f qt=%.3f q3c=%.3f q3s=%.3f qpm=%.3f qcm=%.3f qtm=%.3f q3cm=%.3f q3sm=%.3f\n", cellx, celly, cell->bias, cell->sky, cell->rms, nstar[k], fwmed[k], fwmacro, m2med[k], - qpmed[k], qcmed[k], qtmed[k], - qpmacro, qcmacro, qtmacro); + qpmed[k], qcmed[k], qtmed[k], q3cmed[k], q3smed[k], + qpmacro, qcmacro, qtmacro, q3cmacro, q3smacro); } }