Changeset 23924 for trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c

Timestamp:

Apr 18, 2009, 11:11:15 AM (17 years ago)

Author:

eugene

Message:

import updates from camera group (fits table outputs, fixes to persistent trails); modify import to avoid libpsf (fortran)

File:

• trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (modified) (23 diffs)

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c

@@ -20 +20 @@
 #define EXTERN  /* Define EXTERN to declare variables in params.h */
 #include "burnparams.h"
+#include "persist_fits.h"
 
 int
@@ -28 +29 @@
    HeaderUnit ehu;
    const char* ifilename = "-";
-   const char* ofilename = NULL;
    char extname[FH_MAX_STRLEN+1], otaposn[FH_MAX_STRLEN+1];
    int cellmask[MAXCELL];
@@ -34 +34 @@
    int nextend, cellxy, cell, cellcode;
    int ext, update, restore, psfsize, psfavg;
-   int outfd;
    IMTYPE *buf;
-   const char *burnfile=NULL,  *persistfile=NULL;
+   const char *burnfile=NULL,  *persistfile=NULL, *persistfitsfile=NULL;
+   const char *deltablefitsfile=NULL;
    const char *psffile=NULL, *psfstatfile=NULL;
    CELL OTA[MAXCELL];   /* Cell structure for entire OTA */
 
-   if (argc == 1)  {
-     syntax(argv[0]);
-     exit(EXIT_SUCCESS);
-   }
-
    if(argc > 1 && strncmp(argv[1], "help", 4) == 0) {
-      syntax(argv[0]);
-      exit(EXIT_SUCCESS);
-   }
-
-   if(argc > 1 && strncmp(argv[1], "--help", 6) == 0) {
-      syntax(argv[0]);
-      exit(EXIT_SUCCESS);
-   }
-
-   if(argc > 1 && strncmp(argv[1], "-h", 6) == 0) {
       syntax(argv[0]);
       exit(EXIT_SUCCESS);
@@ -71 +56 @@
       fprintf(stderr, "\rerror: `%s' is not a multi-extension FITS\n",
               ifilename);
+#ifndef JT2DHACK
       exit(EXIT_FAILURE);
+#endif
    }
 
@@ -91 +78 @@
    MAX_READ_NOISE = 20;         /* Maximum believable read noise (ADU) */
    MIN_EADU = 0.3;              /* Minimum believable e/ADU */
+   SAT4SURE  = 60000;           /* Ignore pixels above for noise estimate */
+   MIN_BLAST_PASS = 0.1;        /* Allow blasted cells if they have BIG satfrac */
+   MAX_BLAST_PASS = 0.9;        /* But not if it's all wiped out! */
 
    BURN_THRESH  = 30000;        /* Threshold for onset of burning */
@@ -108 +98 @@
 
    NEGLIGIBLE_TRAIL = 0.5;      /* Don't sweat less than this * sigma */
-
+   EXPIRE_TRAIL_TIME = 2000;    /* Expire a persist after this [sec] */
 
 /* Parse the args */
@@ -155 +145 @@
          burnfile = argv[i] + 9;
 
-/* Output file for burn streaks */
+/* Output text file for burn streaks */
       } else if(strncmp(argv[i], "out=", 4) == 0) {     /* out=fname */
          burnfile = argv[i] + 4;
 
-# if (0)
-/* XXX disable for now: is not yet working */
-/* Alternate Output file for deburned image */
-      } else if(strncmp(argv[i], "outimage=", 9) == 0) {        /* outimage=fname */
-        ofilename = argv[i] + 9;
-# endif
-
-/* Input file for previous burn persistence streaks */
-      } else if(strncmp(argv[i], "trailin=", 8) == 0) { /* trailin=fname */
+/* Input text file for previous burn persistence streaks */
+      } else if(strncmp(argv[i], "infits=", 8) == 0) {  /* in=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 */
+         persistfitsfile = argv[i] + 12;
 
 /* Input file for previous burn persistence streaks */
       } else if(strncmp(argv[i], "in=", 3) == 0) {      /* in=fname */
          persistfile = argv[i] + 3;
+
+/* 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;
 
 /* Output file for PSF gallery */
@@ -244 +235 @@
       } else if(strncmp(argv[i], "thrpsf=", 7) == 0) {  /* thrpsf=thresh */
          if(sscanf(argv[i]+7, "%d", &PSF_THRESH) != 1) {
+            fprintf(stderr, "\rerror: cannot get value from `%s'\n", argv[i]);
+            exit(EXIT_FAILURE);
+         }
+
+      } else if(strncmp(argv[i], "expire=", 7) == 0) {  /* expire=nsec */
+         if(sscanf(argv[i]+7, "%d", &EXPIRE_TRAIL_TIME) != 1) {
             fprintf(stderr, "\rerror: cannot get value from `%s'\n", argv[i]);
             exit(EXIT_FAILURE);
@@ -295 +292 @@
          }
 
+      } else if(strncmp(argv[i], "deltables=", 8) == 0) { /* trailin=fname */
+         deltablefitsfile = argv[i] + 10;        
+
       } else if(strncmp(argv[i], "help", 4) == 0) {     /* help output */
          syntax(argv[0]);
@@ -306 +306 @@
    }
 
+   /* If we're told to remove the table from a FITS file, then
+    * do nothing else. */
+   if(deltablefitsfile) {
+     fprintf(stderr, 
+             "\rRemoving burn tables from %s, all other options ignored.\n",
+             ifilename);
+     if(persist_fits_remove_tables(ihu, deltablefitsfile) != FH_SUCCESS) {
+       exit(EXIT_FAILURE);
+     }
+     else {
+       exit(EXIT_SUCCESS);
+     }
+   }
+   
+   /* If there is no other persistence info supplied, try getting
+    * it from the input FITS file. */
    if(restore && persistfile == NULL) {
-      fprintf(stderr, "\rerror: must specify an input file for restore\n");
-      exit(EXIT_FAILURE);
+     if(persistfitsfile == NULL) persistfitsfile = ifilename;
    }
 
@@ -322 +337 @@
       if(persist_read(OTA, persistfile)) exit(EXIT_FAILURE);
    }
+   else if(persistfitsfile != NULL) {
+     if(persist_fits_read(OTA, persistfitsfile) != FH_SUCCESS) exit(EXIT_FAILURE);
+   }
 
 /* Which OTA is this??? */
@@ -333 +351 @@
    }
 
-   if (ofilename) {
-     outfd = creat (ofilename, S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH);
-     if (outfd == -1) {
-         fprintf(stderr, "\rerror: Cannot open %s for output\n", ofilename);
-         exit(EXIT_FAILURE);
-      }
-
-     if (fh_write(ihu, outfd)) {
-         fprintf(stderr, "\rerror: Trouble writing PHU to %s\n", ofilename);
-         exit(EXIT_FAILURE);
-      }
+   /* Warn that what you're doing might not be a good idea. */
+   {
+     fh_bool burn_applied;
+     
+     if(fh_get_bool(ihu, PHU_NAME_BURN_APPLIED, 
+                    &burn_applied) != FH_SUCCESS) {
+        if(VERBOSE > 0) {
+           fprintf(stderr, 
+                   "warning: Unable to determine whether burn correction "
+                   "already applied - unable to find %s in primary header\n",
+                   PHU_NAME_BURN_APPLIED);
+        }
+        burn_applied = FH_FALSE;
+     }
+     if(restore && (burn_applied == FH_FALSE)) {
+       fprintf(stderr, 
+               "warning: Restoring old burns, but header indicates no burns previously corrected.\n");             
+     }
+     else if (update && (burn_applied == FH_TRUE)) {
+       fprintf(stderr, 
+               "warning: Applying burn correction, but header indicates burns previously corrected.\n");             
+       }
    }
 
 /* Look at all the MEF's extensions */
+#ifndef JT2DHACK
    for (ext = 1; ext <= nextend; ext++) {
+#else
+   for (ext = 1; ext <= MAX(1,nextend); ext++) {
+#endif
       int naxis, naxis1, naxis2, naxis3;
       int prescan1, ovrscan1, ovrscan2, pontime;
@@ -353 +386 @@
       char xtension[FH_MAX_STRLEN + 1];
 
+#ifdef JT2DHACK
+      if(nextend == 0) {
+         ehu = ihu;
+         sprintf(extname, "%s", "xy00");
+      } else {
+#endif
       if (!(ehu = fh_ehu(ihu, ext)) || 
           fh_get_str(ehu, "EXTNAME", extname, sizeof(extname)) != FH_SUCCESS) {
@@ -377 +416 @@
          fprintf(stderr,
                  "warning: Skipping non-image extension %s\n", extname);
+#ifdef JT2DHACK
+      }
+#endif
 /*
  * %%% When reading from a pipe, does the data still need
@@ -448 +490 @@
       if (fh_get_BZERO(ehu, &bzero_d) == FH_SUCCESS) BZERO = bzero_d;
 
+      if(VERBOSE & VERB_NORM) {
+         printf("nx=%d ny=%d prex=%d postx=%d posty=%d BZERO=%.1f\n", 
+                naxis1, naxis2, prescan1, ovrscan1, ovrscan2, BZERO);
+      }
+
       buf = (IMTYPE*)malloc(naxis1*naxis2*naxis3*sizeof(short));
       if (fh_read_padded_image(ehu, fh_file_desc(ehu), buf,
@@ -471 +518 @@
 /* Get bias and sky levels */
          err = cell_stats(nx, ny, naxis1, naxis2, imbuf, OTA+cell);
+         if(err) {
+            if(VERBOSE > 0) {
+               fprintf(stderr, "logonly: error getting bias/sky/rms for cell %d, skipping...\n", cell);
+            }
+            continue;
+         }
+
+/* Does this cell get a pass because it's heavily blasted? */
+/* Turn it in to a relaxation factor for the "read noise" */
+         i = OTA[cell].satfrac > MIN_BLAST_PASS &&
+             OTA[cell].satfrac < MAX_BLAST_PASS &&
+             OTA[cell].sky < TRAIL_THRESH ? 30 : 1;
 
 /* Does this cell look kosher? */
          if(OTA[cell].rms*OTA[cell].rms > 
-            OTA[cell].sky/MIN_EADU + MAX_READ_NOISE*MAX_READ_NOISE) {
+            OTA[cell].sky/MIN_EADU + i*i*MAX_READ_NOISE*MAX_READ_NOISE) {
             if(VERBOSE > 0) {
                fprintf(stderr, "logonly: cell %d is unreasonably noisy, skipping...\n", cell);
+               fprintf(stderr, "logonly: bias = %d sky = %d rms = %d satfrac = %.2f\n", 
+                       OTA[cell].bias, OTA[cell].sky, OTA[cell].rms, OTA[cell].satfrac);
             }
             continue;
@@ -506 +567 @@
 
          } else {
-
 /* Restore the old burns */
             burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, 
@@ -514 +574 @@
 /* Write the corrected data back to the FITS. */
          if(update) {
-           if (ofilename) {
-             fh_ehu(ehu, 0);    /* Seek back to the start of data */
-
-             if (fh_write(ehu, outfd)) {
-               fprintf(stderr, "\rerror: Trouble writing EXT %d to %s\n", ext, ofilename);
-               exit(EXIT_FAILURE);
-             }
-
-             if (fh_write_padded_image(ehu, outfd, buf,
-                                       naxis1*naxis2*naxis3*sizeof(short),
-                                       FH_TYPESIZE_16) != FH_SUCCESS) {
+            fh_ehu(ehu, 0);     /* Seek back to the start of data */
+            if (fh_write_padded_image(ehu, fh_file_desc(ehu), buf,
+                                      naxis1*naxis2*naxis3*sizeof(short),
+                                      FH_TYPESIZE_16) != FH_SUCCESS) {
                fprintf(stderr, "\rerror: failed to re-write image data for extension `%s'.\n",
                        extname);
                free(buf);
                exit(EXIT_FAILURE);
-             }
-           } else {
-             fh_ehu(ehu, 0);    /* Seek back to the start of data */
-             if (fh_write_padded_image(ehu, fh_file_desc(ehu), buf,
-                                       naxis1*naxis2*naxis3*sizeof(short),
-                                       FH_TYPESIZE_16) != FH_SUCCESS) {
-               fprintf(stderr, "\rerror: failed to re-write image data for extension `%s'.\n",
-                       extname);
-               free(buf);
-               exit(EXIT_FAILURE);
-             }
-           }
-         }
-      }
-      
+            }
+         }
+      }
+
       free(buf);
    }
 
-   if (ofilename) {
-     if (close(outfd)) {
-       fprintf(stderr, "\rerror: trouble writing data to output file %s\n", ofilename);
-       exit(EXIT_FAILURE);
-     }
-   }
-      
 /* Dump out the postage stamp file */
    if(psffile != NULL) {
@@ -564 +599 @@
    }
 
+/* Write burn info to FITS file. */
+   if(update) persist_fits_write(OTA, ihu);
+
+   
+   if(restore) {    
+     /* 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_destroy(ihu);
 
 /* Write the persistence data for the next image */
    if(burnfile != NULL) persist_write(OTA, burnfile);
+
 
    exit(EXIT_SUCCESS);
@@ -613 +660 @@
                       CELL *cell)
 {
-   int i, j, k, n;
+   int i, j, k, n, nsat;
 
 /* Get bias stats */
@@ -624 +671 @@
 
 /* Get sky stats */
-   for(k=n=0; k<nx*ny; k+=((617*nx)/1000)) {
+   for(k=n=nsat=0; k<nx*ny; k+=((617*nx)/1000)) {
       i = k % nx;
       j = k / nx;
-      median_buf[n++] = data[i+NX*j];
-   }
+      if(data[i+NX*j] > SAT4SURE) {
+         nsat++;
+      } else if(data[i+NX*j] != NODATA) {
+         median_buf[n++] = data[i+NX*j];
+      }
+   }
+
+   if(n < 20) {         /* Better have hit at least 20! */
+      cell->sky = cell->rms = 0.0;
+      return(-1);
+   }
+   cell->satfrac = ((double)nsat) / n;
+
 /* First pass at sky and quartile */
    cell->sky = int_median(n, median_buf);
@@ -656 +714 @@
    printf(" restore={t|f}  Restore the input MEF by adding input fits?\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");
+   printf("                   in table extensions).  If both this and the 'in' input\n");
+   printf("                   file option are specified, then 'in' takes precedence.\n");
    printf(" out=fname      Output file for burn streaks\n");
-   printf(" outimage=imname Output file for MEF image\n");
+
 //   printf(" trailin=fname  Input file for previous burn persistence streaks\n");
+//   printf(" trailinfits=fname  Input FITS file for previous burn persistence streaks\n");
 //   printf(" trailout=fname Output file for burn streaks\n");
+   printf(" deltables=fname Copy mef_file to new FITS file 'fname' with burn streak\n");
+   printf("                   tables removed. NOTE: if specified, all other options\n");
+   printf("                   will be ignored!\n");
    printf(" psf=fname      Output file for PSF FITS stamp gallery\n");
    printf(" psfstat=fname  Output file for PSF statistics listing\n");
@@ -676 +741 @@
    printf(" rmask=X        Diameter growth factor of burned spots\n");
    printf(" bmask=X        Box size growth of burn/star boxes\n");
+   printf(" expire=N       Retire a blasted burn after N seconds\n");
    printf(" quiet={t|f}    Quiet?\n");
    printf(" verbose=N      Set verbosity bits:\n");