Changeset 26766 for trunk/extsrc/gpcsw/gpcsrc/fits/burntool

Timestamp:

Feb 3, 2010, 1:00:09 PM (17 years ago)

Author:

watersc1

Message:

New burntool, with trails retained even if they had a poor slope.

Location:

trunk/extsrc/gpcsw/gpcsrc/fits/burntool

Files:

• burnfix.c (modified) (2 diffs)
• burntool.c (modified) (17 diffs)
• burntool.h (modified) (3 diffs)
• burnutils.c (modified) (1 diff)
• man/burntool.1 (modified) (2 diffs)
• persist_fits.c (modified) (48 diffs)
• persist_fits.h (modified) (1 diff)
• persistfix.c (modified) (1 diff)
• persistio.c (modified) (10 diffs)
• psfstamp.c (modified) (9 diffs)
• stardetect.c (modified) (1 diff)
• trailfit.c (modified) (4 diffs)

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

@@ -109 +109 @@
 /* Restore all the burns */
    for(k=0; k<cell->npersist; k++) {
-      if(!cell->persist[k].fiterr) 
+      if(!cell->persist[k].fiterr && 
+         (cell->persist[k].func == BURN_PWR || 
+          cell->persist[k].func == BURN_EXP) ) {
          sub_fit(nx, ny, NX, buf, &(cell->persist[k]), -1);
+      }
    }
    return(0);
@@ -124 +127 @@
 /* Restore all the burns */
    for(k=0; k<cell->npersist; k++) {
-      if(!cell->persist[k].fiterr) 
+      if(!cell->persist[k].fiterr && 
+         (cell->persist[k].func == BURN_PWR || 
+          cell->persist[k].func == BURN_EXP) ) {
          sub_fit(nx, ny, NX, buf, &(cell->persist[k]), +1);
+      }
    }
    return(0);

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

@@ -50 +50 @@
       fprintf(stderr, "\rerror: burntool could not open file `%s'\n",
               ifilename);
-      exit(EXIT_FAILURE);
+      exit(-314);
    }
    nextend = fh_extensions(ihu);
@@ -57 +57 @@
               ifilename);
 #ifndef JT2DHACK
-      exit(EXIT_FAILURE);
+      exit(-315);
 #endif
    }
@@ -343 +343 @@
              ifilename);
      if(persist_fits_remove_tables(ihu, deltablefitsfile) != FH_SUCCESS) {
-       exit(EXIT_FAILURE);
-     }
-     else {
+       exit(-316);
+     } else {
        exit(EXIT_SUCCESS);
      }
@@ -365 +364 @@
 /* Read the persistence data for this OTA */
    if(persistfile != NULL) {                    /* Text data file */
-      if(persist_read(OTA, persistfile)) exit(EXIT_FAILURE);
+      if(persist_read(OTA, persistfile, apply)) exit(-317);
    } else if(persistfitsfile != NULL) {         /* FITS table */
-      if(persist_fits_read(OTA, persistfitsfile) != FH_SUCCESS)
-         exit(EXIT_FAILURE);
+      if(persist_fits_read(OTA, persistfitsfile, apply) != FH_SUCCESS)
+         exit(-318);
    }
 
@@ -427 +426 @@
                  "\rerror: Cannot read EXTNAME from `%s' for extension #%d\n",
                  ifilename, ext);
-         exit(EXIT_FAILURE);
+         exit(-319);
       }
       if(fh_get_str(ehu, "XTENSION", xtension, sizeof(xtension)) != FH_SUCCESS) {
@@ -433 +432 @@
                  "\rerror: Cannot read XTENSION from `%s' for extension #%d\n",
                  ifilename, ext);
-         exit(EXIT_FAILURE);     
+         exit(-320);
       }
       if(!strcmp("TABLE", xtension)) {
@@ -461 +460 @@
          fprintf(stderr,
                  "\rerror: Illegal cell number %d from '%s'\n", cell, extname);
-         exit(EXIT_FAILURE);
+         exit(-321);
       }
 
@@ -487 +486 @@
           (naxis >= 3 && fh_get_NAXIS3(ehu, &naxis3) != FH_SUCCESS)) {
          fprintf(stderr, "\rerror: Cannot get NAXIS*'\n");
-         exit(EXIT_FAILURE);
+         exit(-322);
       }
       if (naxis != 2) {
@@ -503 +502 @@
          }
          fprintf(stderr, "\rerror: 32bpp support not yet implemented in burntool.\n");
-         exit(EXIT_FAILURE);
+         exit(-324);
       }
 /* Check BSCALE and warn if it is anything other than 1.0.
@@ -525 +524 @@
       }
 
-      buf = (IMTYPE*)malloc(naxis1*naxis2*naxis3*sizeof(short));
+      if( (buf = (IMTYPE*)malloc(naxis1*naxis2*naxis3*sizeof(short))) == NULL) {
+         fprintf(stderr, "\rerror: failed to alloc FITS buffer\n");
+         exit(-325);
+      }
       if (fh_read_padded_image(ehu, fh_file_desc(ehu), buf,
                                naxis1*naxis2*naxis3*sizeof(short),
@@ -531 +533 @@
          fprintf(stderr, "\rerror: failed to read image data for extension `%s'.\n",
                  extname);
-         free(buf);
-         exit(EXIT_FAILURE);
+//       free(buf);
+         exit(-326);
       }
 
@@ -617 +619 @@
                fprintf(stderr, "\rerror: failed to re-write image data for extension `%s'.\n",
                        extname);
-               free(buf);
-               exit(EXIT_FAILURE);
+//             free(buf);
+               exit(-327);
             }
          }
@@ -667 +669 @@
    if(2*SKY_MARG*ny > nmedian_buf) {
       if(median_buf != NULL) free(median_buf);
-      median_buf = (int *)calloc(2*SKY_MARG*ny, sizeof(int));
+      if( (median_buf = (int *)calloc(2*SKY_MARG*ny, sizeof(int))) == NULL) {
+         fprintf(stderr, "\rerror: failed to alloc median buffer\n");
+         exit(-667);
+      }
+
       nmedian_buf = 2*SKY_MARG*ny;
    }
@@ -674 +680 @@
    if(NX*NY > nimbuf) {
       if(imbuf != NULL) free(imbuf);
-      imbuf = (int *)calloc(NX*NY, sizeof(DTYPE));
+      if( (imbuf = (int *)calloc(NX*NY, sizeof(DTYPE))) == NULL) {
+         fprintf(stderr, "\rerror: failed to alloc cell copy\n");
+         exit(-668);
+      }
       nimbuf = NX*NY;
    }
@@ -681 +690 @@
    if(NX*NY > nmbuf) {
       if(mbuf != NULL) free(mbuf);
-      mbuf = (int *)calloc(NX*NY, sizeof(MTYPE));
+      if( (mbuf = (int *)calloc(NX*NY, sizeof(MTYPE))) == NULL) {
+         fprintf(stderr, "\rerror: failed to alloc mask buffer\n");
+         exit(-669);
+      }
       nmbuf = NX*NY;
    }
@@ -688 +700 @@
    if(NX*NY > nmsbuf) {
       if(msbuf != NULL) free(msbuf);
-      msbuf = (int *)calloc(NX*NY, sizeof(MTYPE));
+      if( (msbuf = (int *)calloc(NX*NY, sizeof(MTYPE))) == NULL) {
+         fprintf(stderr, "\rerror: failed to alloc veto buffer\n");
+         exit(-670);
+      }
       nmsbuf = NX*NY;
    }
@@ -755 +770 @@
    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(" tableonly={t|f} Calculate fits but do *not* modify the input MEF pixels, 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");

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

@@ -91 +91 @@
 #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 */
 
@@ -175 +176 @@
 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_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,
@@ -208 +213 @@
 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);

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

@@ -290 +290 @@
    return(0);
 }
+
+/****************************************************************/
+/* Return of RMS relative to Fit y = ax + b */
+STATIC int linearrms(int npt, double *x, double *y, double a, double b, 
+                     double *rms)
+{
+   *rms = 0.0;
+   while(--npt >= 0) *rms += (y[npt]-a*x[npt]-b)*(y[npt]-a*x[npt]-b);
+   if(*rms > 0) *rms = sqrt(*rms);
+   return(0);
+}

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1

@@ -60 +60 @@
         correction is applied) and they will propagate from input to output
         until they finally achieve a legal fit with negligible amplitude.
+        Note that if burntool decides that a fit had a bad slope but was 
+        non-negligible it writes a function type "BURN_POSSLOPE" = 4 so
+        the box is available for masking if desired, and can be refitted
+        on subsequent images.
 
         Burntool also identifies really blasted areas which are saturated from
@@ -227 +231 @@
 
         tableonly={t|f}
-                Calculate fits but do *not* modify the input MEF images, 
+                Calculate fits but do *not* modify the input MEF pixels, 
                 only write the fit data as a FITS and/or text table (default f)
 

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

@@ -52 +52 @@
    AREA_TABLE_COL_SXFIT,
    AREA_TABLE_COL_EXFIT,
+   AREA_TABLE_COL_FITERR,
    
    /* Add new columns above this line. */
@@ -86 +87 @@
    { "sxfit", "Starting column for fit",              "pixels",  FH_TABLE_FORMAT_INT,    3,     0 },
    { "exfit", "Ending column for fit",                "pixels",  FH_TABLE_FORMAT_INT,    3,     0 },
+   { "fiterr", "Error code of fit",                "",  FH_TABLE_FORMAT_INT,    3,     0 },
 };
 
@@ -222 +224 @@
       for(k=0; k<cell[j].npersist; k++) 
       {
+/* Retire old burns */
+         if(cell[j].time - cell[j].persist[k].time > EXPIRE_TRAIL_TIME)
+            continue;
+
          if(PERSIST_RETAIN) {
 /* Keep fits which have a dubious slope */
@@ -306 +312 @@
       (fh_table_write_value(table, data, row, AREA_TABLE_COL_NFIT,  &(area->nfit))  != FH_SUCCESS) ||
       (fh_table_write_value(table, data, row, AREA_TABLE_COL_SXFIT, &(area->sxfit)) != FH_SUCCESS) ||
-      (fh_table_write_value(table, data, row, AREA_TABLE_COL_EXFIT, &(area->exfit)) != FH_SUCCESS))
-   {
-      fprintf(stderr, "error: Error writing data to row %d of area table.\n", row);
+      (fh_table_write_value(table, data, row, AREA_TABLE_COL_EXFIT, &(area->exfit)) != FH_SUCCESS) ||
+      (fh_table_write_value(table, data, row, AREA_TABLE_COL_FITERR, &(area->fiterr)) != FH_SUCCESS))
+   {
+      fprintf(stderr, "\rerror: Error writing data to row %d of area table.\n", row);
       return FH_BAD_VALUE;
    }
@@ -346 +353 @@
       for(k=0; k<cell[j].npersist; k++) 
       {
+/* Retire old burns */
+         if(cell[j].time - cell[j].persist[k].time > EXPIRE_TRAIL_TIME)
+            continue;
+
          if(PERSIST_RETAIN) {
 /* Keep fits which have a dubious slope */
@@ -422 +433 @@
       (fh_table_write_value(table, data, row, FIT_TABLE_COL_ZERO, &zero) != FH_SUCCESS))
    {
-      fprintf(stderr, "error: Error writing data to row %d of fit table.\n", row);
+      fprintf(stderr, "\rerror: Error writing data to row %d of fit table.\n", row);
       return FH_BAD_VALUE;
    }
@@ -460 +471 @@
       for(k=0; k<cell[j].npersist; k++) 
       {
+/* Retire old burns */
+         if(cell[j].time - cell[j].persist[k].time > EXPIRE_TRAIL_TIME)
+            continue;
+
          if(PERSIST_RETAIN) {
 /* Keep fits which have a dubious slope */
@@ -526 +541 @@
  */
 static fh_result
-read_area(fhTable * table, void * data, int row)
+read_area(fhTable * table, void * data, int row, int apply)
 {
    int cell_num;
@@ -534 +549 @@
    {
       fprintf(stderr,
-              "error: Unable to get cell number from row %d of burn area table\n",
+              "\rerror: Unable to get cell number from row %d of burn area table\n",
               row);
       return FH_BAD_VALUE;
@@ -542 +557 @@
    if((cell_num < 0) || (cell_num > MAXCELL))
    {
-      fprintf(stderr, "error: illegal cell %d in area table row %d\n",
+      fprintf(stderr, "\rerror: illegal cell %d in area table row %d\n",
               cell_num, row);
       boxbuf[row].cell = -1;
@@ -573 +588 @@
       (fh_table_read_value(table, data, row, AREA_TABLE_COL_NFIT,  &(boxbuf[row].nfit)) != FH_SUCCESS) ||
       (fh_table_read_value(table, data, row, AREA_TABLE_COL_SXFIT, &(boxbuf[row].sxfit)) != FH_SUCCESS) ||
-      (fh_table_read_value(table, data, row, AREA_TABLE_COL_EXFIT, &(boxbuf[row].exfit)) != FH_SUCCESS))
+      (fh_table_read_value(table, data, row, AREA_TABLE_COL_EXFIT, &(boxbuf[row].exfit)) != FH_SUCCESS) ||
+      (fh_table_read_value(table, data, row, AREA_TABLE_COL_FITERR, &(boxbuf[row].fiterr)) != FH_SUCCESS))
    {
       fprintf(stderr,
-              "error: Error reading values from row %d of burn area table\n",
+              "\rerror: Error reading values from row %d of burn area table\n",
               row);
       boxbuf[row].nfit = 0;
@@ -593 +609 @@
       boxbuf[row].xfit = (int *)calloc(boxbuf[row].nfit, sizeof(int));
       boxbuf[row].yfit = (int *)calloc(boxbuf[row].nfit, sizeof(int));
-      boxbuf[row].fiterr = 0;
+
+      if(boxbuf[row].zero == NULL ||
+         boxbuf[row].xfit == NULL ||
+         boxbuf[row].yfit == NULL) {
+         fprintf(stderr, "\rerror: failed to alloc boxbuf\n");
+         exit(-671);
+      }
+// 100203 JT: fiterr now saved and read, refit if not just an "apply"
+      if(!apply) boxbuf[row].fiterr = 0;
    }
 
@@ -616 +640 @@
  */
 static fh_result
-read_area_table(HeaderUnit hu, fhTable * table)
+read_area_table(HeaderUnit hu, fhTable * table, int apply)
 {
    fh_result result = FH_INVALID;
@@ -627 +651 @@
       (fh_get_int(hu, "TFIELDS", &num_cols) != FH_SUCCESS))
    {
-      fprintf(stderr, "error: Unable to find required keywords for area table dimensions\n");
+      fprintf(stderr, "\rerror: Unable to find required keywords for area table dimensions\n");
       return FH_NOT_FOUND;
    }
@@ -634 +658 @@
    if(num_cols != table->num_cols) 
    {
-      fprintf(stderr, "error: %d-column area table found, expected %d cols.\n",
+      fprintf(stderr, "\rerror: %d-column area table found, expected %d cols.\n",
               num_cols, table->num_cols);
       return FH_BAD_VALUE;
@@ -642 +666 @@
    {
       fprintf(stderr, 
-              "error: too many boxes in area table. Max is %d, got %d\n",
+              "\rerror: too many boxes in area table. Max is %d, got %d\n",
               MAXBURN, num_rows);
       return FH_BAD_VALUE;
@@ -654 +678 @@
    {
       fprintf(stderr, 
-              "error: Unable to allocate %d bytes for area table.\n",
+              "\rerror: Unable to allocate %d bytes for area table.\n",
               table->table_size);
       return FH_NO_MEMORY;
@@ -662 +686 @@
                     table->table_size, FH_TYPESIZE_8) != FH_SUCCESS)
    {
-      fprintf(stderr, "error: Unable to map area table body for reading.\n");
+      fprintf(stderr, "\rerror: Unable to map area table body for reading.\n");
       free(data);
       return FH_INVALID;
@@ -669 +693 @@
    for(i = 0; i < num_rows; i++)
    {
-      if((result = read_area(table, data, i)) != FH_SUCCESS) break;
+      if((result = read_area(table, data, i, apply)) != FH_SUCCESS) break;
    }
 
@@ -713 +737 @@
    {
       fprintf(stderr,
-              "error: Unable to get cell number from row %d of fit table\n",
+              "\rerror: Unable to get cell number from row %d of fit table\n",
               row);
       return FH_BAD_VALUE;
@@ -720 +744 @@
    if((*cell < 0) || (*cell > MAXCELL))
    {
-      fprintf(stderr, "error: illegal cell %d in fit table row %d\n",
+      fprintf(stderr, "\rerror: illegal cell %d in fit table row %d\n",
               *cell, row);
       
@@ -733 +757 @@
    {
       fprintf(stderr,
-              "error: Error reading values from row %d of burn area table\n",
+              "\rerror: Error reading values from row %d of burn area table\n",
               row);
       return FH_BAD_VALUE;
@@ -782 +806 @@
       (fh_get_int(hu, "TFIELDS", &num_cols) != FH_SUCCESS))
    {
-      fprintf(stderr, "error: Unable to find required keywords for fit table dimensions\n");
+      fprintf(stderr, "\rerror: Unable to find required keywords for fit table dimensions\n");
       return FH_NOT_FOUND;
    }
@@ -789 +813 @@
    if(num_cols != table->num_cols) 
    {
-      fprintf(stderr, "error: %d-column fit table found, expected %d cols.\n",
+      fprintf(stderr, "\rerror: %d-column fit table found, expected %d cols.\n",
               num_cols, table->num_cols);
       return FH_BAD_VALUE;
@@ -805 +829 @@
    {
       fprintf(stderr, 
-              "error: Unable to allocate %d bytes for area table.\n",
+              "\rerror: Unable to allocate %d bytes for area table.\n",
               table->table_size);
       return FH_NO_MEMORY;
@@ -811 +835 @@
    if(fh_read_image(hu, fh_file_desc(hu), data, table->table_size, FH_TYPESIZE_8) != FH_SUCCESS)
    {
-      fprintf(stderr, "error: Unable to read fit table body.\n");
+      fprintf(stderr, "\rerror: Unable to read fit table body.\n");
       free(data);
       return FH_INVALID;
@@ -836 +860 @@
          {
             fprintf(stderr, 
-                    "error: Fit in table row %d does not match current area "
+                    "\rerror: Fit in table row %d does not match current area "
                     "(table cell=%d cx=%d cy=%d vs. area cell=%d cx=%d cy=%d)\n",
                     fit_table_row,
@@ -882 +906 @@
  */
 fh_result
-persist_fits_read(CELL *cell, const char * filename)
+persist_fits_read(CELL *cell, const char * filename, int apply)
 {
    HeaderUnit phu;
@@ -891 +915 @@
    {
       fprintf(stderr, 
-              "error: unable to create header structure for FITS file \"%s\"\n",
+              "\rerror: unable to create header structure for FITS file \"%s\"\n",
               filename);
       return FH_NO_MEMORY;
@@ -899 +923 @@
    if ((result = fh_file(phu, filename, FH_FILE_RDONLY)) != FH_SUCCESS)
    {
-      fprintf(stderr, 
-              "error: unable to open FITS file \"%s\"\n",
+      fprintf(stderr, "\rerror: unable to open FITS file \"%s\"\n",
               filename);
       fh_destroy(phu);
@@ -911 +934 @@
    if(!area_hu || !fit_hu) 
    {
-      fprintf(stderr, 
-              "error: Unable to find persistence info in FITS file.\n");
+      fprintf(stderr, "\rerror: Unable to find persistence info in FITS file.\n");
       fh_destroy(phu);
       return FH_NOT_FOUND;
@@ -921 +943 @@
 
    fh_ehu_by_extname(phu, DEFAULT_EXTNAME_AREA_TABLE);
-   if((result = read_area_table(area_hu, &area_table)) != FH_SUCCESS)
+   if((result = read_area_table(area_hu, &area_table, apply)) != FH_SUCCESS)
    {
       free(area_table.strbuf);
       free(fit_table.strbuf);
       fh_destroy(phu);
+      fprintf(stderr, "\rerror: Unable to read area table from FITS file.\n");
       return result;
    }
@@ -935 +958 @@
       free(fit_table.strbuf);
       fh_destroy(phu);
+      fprintf(stderr, "\rerror: Unable to read fit table from FITS file.\n");
       return result;
    }
@@ -961 +985 @@
       for(k=0; k<MAXCELL; k++) {
          if( (i=cell[k].npersist) > 0) {
-            cell[k].persist = (OBJBOX *)calloc(i, sizeof(OBJBOX));
+            if( (cell[k].persist = (OBJBOX *)calloc(i, sizeof(OBJBOX))) == NULL) {
+               fprintf(stderr, "\rerror: failed to alloc cell persist buffer\n");
+               exit(-672);
+            }
             cell[k].npersist = 0;
          }
@@ -1018 +1045 @@
    {
       fprintf(stderr, 
-              "error: Unable to write correction info to FITS file. Correction FITS tables already exist.\n");
+              "\rerror: Unable to write correction info to FITS file. Correction FITS tables already exist.\n");
       return FH_NO_SPACE;
    }
@@ -1025 +1052 @@
    if ((fd = fh_file_desc(phu)) == -1)
    {
-      fprintf(stderr, "error: header passed has no associated file.\n");
+      fprintf(stderr, "\rerror: header passed has no associated file.\n");
       return FH_INVALID;
    }
@@ -1033 +1060 @@
    {
       fprintf(stderr,
-              "error: Unable to seek to end of file.\n");
+              "\rerror: Unable to seek to end of file.\n");
       return FH_IN_ERRNO;
    }
@@ -1042 +1069 @@
    {
       fprintf(stderr, 
-              "error: Unable to create headers for correction FITS tables.\n");
+              "\rerror: Unable to create headers for correction FITS tables.\n");
       if(area_hu) fh_destroy(area_hu);
       if(fit_hu) fh_destroy(fit_hu);
@@ -1066 +1093 @@
       (fh_munmap_table(area_hu) != FH_SUCCESS))
    {
-      fprintf(stderr, "error: Error encountered writing area table to FITS file.\n");
+      fprintf(stderr, "\rerror: Error encountered writing area table to FITS file.\n");
       fh_destroy(fit_hu);
       free(area_table.strbuf);
@@ -1087 +1114 @@
       (fh_munmap_table(fit_hu) != FH_SUCCESS))
    {
-      fprintf(stderr, "error: Error encountered writing area table to FITS file.\n");
+      fprintf(stderr, "\rerror: Error encountered writing area table to FITS file.\n");
       free(area_table.strbuf);
       free(fit_table.strbuf);
@@ -1149 +1176 @@
    if(fh_get_int(phu_in, "NEXTEND", &num_extensions) != FH_SUCCESS)
    {
-      fprintf(stderr, "error: Unable to get NEXTEND from input primary FITS header.\n");
+      fprintf(stderr, "\rerror: Unable to get NEXTEND from input primary FITS header.\n");
       return FH_INVALID;
    }
@@ -1157 +1184 @@
    {
       fprintf(stderr, 
-              "error: Unable to find persistence info in FITS file.\n");
+              "\rerror: Unable to find persistence info in FITS file.\n");
       return FH_INVALID;
    }
@@ -1175 +1202 @@
    if((fd_out = open(fileout, O_CREAT | O_RDWR, 0644)) < 0)
    {
-      fprintf(stderr, "error: Failed to open \"%s\" for output.\n",
+      fprintf(stderr, "\rerror: Failed to open \"%s\" for output.\n",
               fileout);
       exit(EXIT_FAILURE);
@@ -1192 +1219 @@
       (fh_write(phu_out, fd_out) != FH_SUCCESS))
    {
-      fprintf(stderr, "error: Unable to copy primary header to \"%s\"\n",
+      fprintf(stderr, "\rerror: Unable to copy primary header to \"%s\"\n",
               fileout);
       fh_destroy(phu_out);
@@ -1209 +1236 @@
       if(!(ehu_in = fh_ehu(phu_in, i)))
       {
-         fprintf(stderr, "error: Unable to read extension %d from input FITS file\n", i);
+         fprintf(stderr, "\rerror: Unable to read extension %d from input FITS file\n", i);
          fh_destroy(phu_out);
          close(fd_out);
@@ -1217 +1244 @@
       if (fh_get_str(ehu_in, "EXTNAME", extname, sizeof(extname)) != FH_SUCCESS)
       {
-         fprintf(stderr, "error: Unable to get EXTNAME from extension %d in input FITS file.\n", i); 
+         fprintf(stderr, "\rerror: Unable to get EXTNAME from extension %d in input FITS file.\n", i); 
          fh_destroy(phu_out);
          close(fd_out);
@@ -1233 +1260 @@
             (fh_reserve(ehu_out, (2880/80)) != FH_SUCCESS))
          {
-            fprintf(stderr, "error: Unable to copy extension %d to %s.\n", i, fileout); 
+            fprintf(stderr, "\rerror: Unable to copy extension %d to %s.\n", i, fileout); 
             fh_destroy(phu_out);
             close(fd_out);
@@ -1248 +1275 @@
             (fh_copy_padded_image(ehu_out, fd_out, fh_file_desc(ehu_in)) != FH_SUCCESS))
          {
-            fprintf(stderr, "error: Unable to copy extension %d to %s.\n", i, fileout); 
+            fprintf(stderr, "\rerror: Unable to copy extension %d to %s.\n", i, fileout); 
             fh_destroy(phu_out);
             close(fd_out);

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persist_fits.h

@@ -24 +24 @@
 
 fh_result
-persist_fits_read(CELL *cell, const char * filename);
+persist_fits_read(CELL *cell, const char * filename, int apply);
 
 fh_result

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

@@ -34 +34 @@
          }
          continue;
+      }
+
+/* This had a significant positive slope, what do we do now? */
+      if( (cell->persist)[k].func == BURN_POSSLOPE) {
+/* try again, so let it slide through... */
       }
 

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

@@ -16 +16 @@
 /****************************************************************/
 /* persist_read(): Read all the persistence trails from a file */
-STATIC int persist_read(CELL *cell, const char *infile)
+STATIC int persist_read(CELL *cell, const char *infile, int apply)
 {
    int i, k, nbox=0;
@@ -38 +38 @@
    while(fgets(line, 1024, fp) != NULL) {
       if(line[0] == '#') continue;
-      sscanf(line, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %lf %d %d %d",
+      sscanf(line, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %lf %d %d %d %d",
              &boxbuf[nbox].cell, &boxbuf[nbox].time,
              &boxbuf[nbox].cx, &boxbuf[nbox].cy,
@@ -50 +50 @@
              &boxbuf[nbox].func, &boxbuf[nbox].up,
              &boxbuf[nbox].slope, &boxbuf[nbox].nfit,
-             &boxbuf[nbox].sxfit, &boxbuf[nbox].exfit);
+             &boxbuf[nbox].sxfit, &boxbuf[nbox].exfit,
+             &boxbuf[nbox].fiterr);
       if(boxbuf[nbox].nfit > 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<boxbuf[nbox].nfit; i++) {
             if(fgets(line, 1024, fp) == NULL) {
@@ -64 +72 @@
          }
       }
-      boxbuf[nbox].fiterr = 0;
+// 100203 JT: fiterr now saved and read, refit if not just an "apply"
+      if(!apply) boxbuf[nbox].fiterr = 0;
 /* Augment counts */
       k = boxbuf[nbox].cell;
@@ -85 +94 @@
    for(k=0; k<MAXCELL; k++) {
       if( (i=cell[k].npersist) > 0) {
-         cell[k].persist = (OBJBOX *)calloc(i, sizeof(OBJBOX));
+         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;
       }
@@ -242 +254 @@
       for(k=0; k<cell[j].npersist; k++) {
 
+/* Retire old burns */
+         if(cell[j].time - cell[j].persist[k].time > EXPIRE_TRAIL_TIME)
+            continue;
+
          if(PERSIST_RETAIN) {
 /* Keep fits which have a dubious slope */
@@ -252 +268 @@
             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\n",
+         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, 
@@ -264 +280 @@
                  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].sxfit, cell[j].persist[k].exfit, 
+                 cell[j].persist[k].fiterr);
          for(i=0; i<cell[j].persist[k].nfit; i++) {
             fprintf(fp, "%3d %3d %8.4f\n", cell[j].persist[k].xfit[i], 
@@ -288 +305 @@
 
          i = (cell[j].burn[k].ex - cell[j].burn[k].sx + 1) / 2;
-         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\n", 
+         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].burn[k].time, 
                  cell[j].burn[k].cx, cell[j].burn[k].cy,
@@ -300 +317 @@
                  cell[j].burn[k].func, cell[j].burn[k].up, 
                  cell[j].burn[k].slope, cell[j].burn[k].nfit,
-                 cell[j].burn[k].sxfit, cell[j].burn[k].exfit);
+                 cell[j].burn[k].sxfit, cell[j].burn[k].exfit, 
+                 cell[j].burn[k].fiterr);
          for(i=0; i<cell[j].burn[k].nfit; i++) {
             fprintf(fp, "%3d %3d %8.4f\n", cell[j].burn[k].xfit[i], 

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

@@ -112 +112 @@
 
 /* Make a postage stamp for it */
-      box[k].stamp = (IMTYPE *)calloc(size*size, sizeof(IMTYPE));
+      if( (box[k].stamp = (IMTYPE *)calloc(size*size, sizeof(IMTYPE))) == NULL) {
+         fprintf(stderr, "\rerror: failed to stamp buffer\n");
+         exit(-675);
+      }
+
       sum = 0;
       for(j=y0; j<=y1; j++) {
@@ -136 +140 @@
    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[5], qt, xfp, yfp, pi=4*atan(1.0);
+   double scale, phi, fwhm[3], q[7], qt, xfp, yfp, pi=4*atan(1.0);
    IMTYPE *median_image;
    CELL *cell;
@@ -171 +175 @@
 
 /* Create a median image and write as the first one */
-   median_image = (IMTYPE *)calloc(nx*ny, sizeof(IMTYPE));
+   if( (median_image = (IMTYPE *)calloc(nx*ny, sizeof(IMTYPE))) == NULL) {
+      fprintf(stderr, "\rerror: failed to alloc median PSF image\n");
+      exit(-676);
+   }
    for(i=0; i<nx*ny; i++) {
       nstar = 0;
@@ -256 +263 @@
    int k, l, cellx, celly, ota_xid, ota_yid, nfwave, nqavg;
    double xota, yota, xfp, yfp, phi;
-   double fwhm[3], q[5], fw[MAXPSFMEDIAN];
+   double fwhm[3], q[7], fw[MAXPSFMEDIAN];
    double m2[MAXPSFMEDIAN], qp[MAXPSFMEDIAN], qc[MAXPSFMEDIAN];
    double qt[MAXPSFMEDIAN], fwavg[MAXPSFMEDIAN];
+   double qpavg[MAXPSFMEDIAN], qcavg[MAXPSFMEDIAN], qtavg[MAXPSFMEDIAN];
    double q3c[MAXPSFMEDIAN], q3s[MAXPSFMEDIAN];
-   double qpavg[MAXPSFMEDIAN], qcavg[MAXPSFMEDIAN], qtavg[MAXPSFMEDIAN];
    double q3cavg[MAXPSFMEDIAN], q3savg[MAXPSFMEDIAN];
+   double q1c[MAXPSFMEDIAN], q1s[MAXPSFMEDIAN];
+   double q1cavg[MAXPSFMEDIAN], q1savg[MAXPSFMEDIAN];
    int nstar[MAXCELL], nfw[MAXCELL];
    double fwmed[MAXCELL], m2med[MAXCELL];
    double qpmed[MAXCELL], qcmed[MAXCELL], qtmed[MAXCELL];
    double q3cmed[MAXCELL], q3smed[MAXCELL];
-   double qpmacro, qcmacro, qtmacro, fwmacro, q3cmacro, q3smacro;
+   double q1cmed[MAXCELL], q1smed[MAXCELL];
+   double qpmacro, qcmacro, qtmacro, fwmacro, q3cmacro, q3smacro, q1cmacro, q1smacro;
    FILE *fp;
    CELL *cell;
@@ -310 +320 @@
                      q3c[nstar[k]] = q[3];
                      q3s[nstar[k]] = q[4];
+                     q1c[nstar[k]] = q[5];
+                     q1s[nstar[k]] = q[6];
 /* Get the position in the focal plane and therefore the qt statistic */
                      psc_cell_to_pixel(cellx, celly, 0.5*PSC_HCELL/PSC_PIXEL, 
@@ -333 +345 @@
                   q3cmed[k] = double_median(nstar[k], q3c);
                   q3smed[k] = double_median(nstar[k], q3s);
+                  q1cmed[k] = double_median(nstar[k], q1c);
+                  q1smed[k] = double_median(nstar[k], q1s);
                } else {
                   m2med[k] = qpmed[k] = qcmed[k] = qtmed[k] = -99.99;
                   q3cmed[k] = q3smed[k] = -99.99;
+                  q1cmed[k] = q1smed[k] = -99.99;
                }
 /* Toss these results into the macrocell median hopper */
@@ -348 +363 @@
                      q3cavg[nqavg] = q3c[l];
                      q3savg[nqavg] = q3s[l];
+                     q1cavg[nqavg] = q1c[l];
+                     q1savg[nqavg] = q1s[l];
                      nqavg++;
                   }
@@ -359 +376 @@
          q3cmacro = double_median(nqavg, q3cavg);
          q3smacro = double_median(nqavg, q3savg);
-         if(nqavg == 0) qpmacro = qcmacro = qtmacro = q3cmacro = q3smacro = -99.99;
+         q1cmacro = double_median(nqavg, q1cavg);
+         q1smacro = double_median(nqavg, q1savg);
+         if(nqavg == 0) qpmacro = qcmacro = qtmacro = q3cmacro = q3smacro = q1cmacro = q1smacro = -99.99;
          fwmacro = double_median(nfwave, fwavg);
 
@@ -370 +389 @@
                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 q3c=%.3f q3s=%.3f qpm=%.3f qcm=%.3f qtm=%.3f q3cm=%.3f q3sm=%.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 q1c=%.3f q1s=%.3f qpm=%.3f qcm=%.3f qtm=%.3f q3cm=%.3f q3sm=%.3f q1cm=%.3f q1sm=%.3f\n", 
                        cellx, celly, cell->bias, cell->sky, cell->rms, 
                        nstar[k], fwmed[k], fwmacro, m2med[k], 
-                       qpmed[k], qcmed[k], qtmed[k], q3cmed[k], q3smed[k],
-                       qpmacro, qcmacro, qtmacro, q3cmacro, q3smacro);
+                       qpmed[k], qcmed[k], qtmed[k], q3cmed[k], q3smed[k], 
+                       q1cmed[k], q1smed[k],
+                       qpmacro, qcmacro, qtmacro, q3cmacro, q3smacro, 
+                       q1cmacro, q1smacro);
             }
          }

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

@@ -193 +193 @@
    cell->burn = (OBJBOX *)calloc(cell->nburn, sizeof(OBJBOX));
    cell->star = (OBJBOX *)calloc(cell->nstar, sizeof(OBJBOX));
+   if( cell->burn == NULL || cell->star == NULL) {
+      fprintf(stderr, "\rerror: failed to alloc burn box\n");
+      exit(-678);
+   }
 
 /* Copy the boxes to the cell info structure */

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

@@ -30 +30 @@
 
    if(box->func != BURN_PWR && box->func != BURN_EXP) {
-      fprintf(stderr, "error: unimplemented fit function %d\n", box->func);
+      fprintf(stderr, "error: unimplemented fit function %d (err %d)\n", 
+              box->func, box->fiterr);
       return(-1);
    }
@@ -137 +138 @@
    box->xfit = (int *)calloc(xe-xs+1, sizeof(int));
    box->yfit = (int *)calloc(xe-xs+1, sizeof(int));
+   if(box->zero == NULL || box->xfit == NULL || box->yfit == NULL) {
+      fprintf(stderr, "\rerror: failed to alloc box memory\n");
+      exit(-679);
+   }
 /* Some defaults */
    box->slope = 0.0;
@@ -229 +234 @@
       box->nfit = 0;
       box->fiterr = FIT_SLOPE_ERROR;
+/* Check whether it's a significant trail (but with pos slope) or just noise */
+      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;
       return(-1);
    }
@@ -268 +278 @@
 /* FIXME: what's a really good criterion for negligible fit? */
 /* 100 pixels up fit is zsum or ~zsum/e */
-      if(zsum > NEGLIGIBLE_TRAIL*rms) {
+//      if(zsum > NEGLIGIBLE_TRAIL*rms) {
+      if(zsum > NEGLIGIBLE_TRAIL*rms ||
+/* 100113: but also evaluate at the end of the fit for stubby trails */
+         zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) {
 
 /* Ascertain the starting point of where the fit is good */