Changeset 24522

Timestamp:

Jun 19, 2009, 4:11:25 PM (17 years ago)

Author:

eugene

Message:

updates to tiltystreak from JT & SI

File:

• trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/tiltystreak.c (modified) (15 diffs)

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

@@ +1 @@
+/* -*- c-file-style: "Ellemtel" -*- */
 /* tiltystreak.c - fit and subtract tilted streaks from a row of cells */
 /* Assumes that biastool has already been run -- zero streakiness near bias */
@@ -52 +53 @@
 static int parabola_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]);
 int linearfit(int npt, double *y, int nx, double *x, double *param);
+static void syntax(const char *prog);
 
 static int nmedian_buf=0;
@@ -59 +61 @@
 int main(int argc, const char* argv[])
 {
-   int i, row, col, ext;
+   int i, j, row, col, ext;
    HeaderUnit ihu = fh_create();
    HeaderUnit ehu = NULL;
@@ -66 +68 @@
    int nextend;
    CELL otarow[NCELL];
-   int naxis, naxis1, naxis2, naxis3;
+   int naxis, naxis1=0, naxis2=0;
+   int n1, n2;
    int prescan1, ovrscan1, ovrscan2;
    double bzero_d;
    int bzero = 32768, firstime = 1, algo=2;
+   int cellmask[NCELL*NCELL];
 
    if(argc < 2) {
+      syntax(argv[0]);
       fprintf(stderr, "error: syntax: tiltystreak Uncompressed_GPC1_MEF.fits\n");
       exit(EXIT_FAILURE);
@@ -77 +82 @@
    ifilename = argv[1];
 
-/* Choose different algorithms? */
-   if(argc > 2 && strncmp(argv[2], "alg=", 4) == 0) {
-      if(strcmp(argv[2]+4, "line") == 0) {
-         algo = 1;
-      } else if(strcmp(argv[2]+4, "parabola") == 0) {
-         algo = 2;
-      } else if(strcmp(argv[2]+4, "bent") == 0) {
-         algo = 3;
+/* Defaults */
+   for(i=0; i<NCELL*NCELL; i++) cellmask[i] = 1;        /* Use all cells */
+   algo = 2;                    /* Default algorithm is a parabola */
+
+/* Process arguments */
+   for(i=2; i<argc; i++) {
+
+/* Veto cells by a 64 digit mask (0/1 for cells 0-63) */
+      if(strncmp(argv[i], "mask=", 5) == 0) {   /* mask=64_digits */
+         if(strlen(argv[i]) < 5+NCELL*NCELL) {
+            fprintf(stderr, "\rerror: cannot get mask from `%s'\n", argv[i]);
+            exit(EXIT_FAILURE);
+         } else {
+            for(j=0; j<NCELL*NCELL; j++) cellmask[j] = argv[i][5+j] != '0';
+         }
+
+      } else if(strncmp(argv[2], "alg=", 4) == 0) {     /* Algorithm */
+         if(strcmp(argv[2]+4, "line") == 0) {
+            algo = 1;
+         } else if(strcmp(argv[2]+4, "parabola") == 0) {
+            algo = 2;
+         } else if(strcmp(argv[2]+4, "bent") == 0) {
+            algo = 3;
+         } else {
+            fprintf(stderr, "error: `%s' is not a valid algorithm\n", argv[2]);
+            exit(EXIT_FAILURE);
+         }
+
+      } else if(strncmp(argv[i], "help", 4) == 0) {     /* help output */
+         syntax(argv[0]);
+         exit(EXIT_SUCCESS);
+
       } else {
-         fprintf(stderr, "error: `%s' is not a valid algorithm\n", argv[2]);
+         fprintf(stderr, "\rerror: unrecognized option `%s'\n", argv[i]);
+         syntax(argv[0]);
          exit(EXIT_FAILURE);
       }
    }
+
 
    if (fh_file(ihu, ifilename, FH_FILE_RDWR) != FH_SUCCESS) {
@@ -139 +170 @@
 
 /* Positioned at the desired cell now */
-         naxis3 = 1;
          if (fh_get_PRESCAN1(ehu, &prescan1) != FH_SUCCESS ||
              fh_get_OVRSCAN1(ehu, &ovrscan1) != FH_SUCCESS ||
              fh_get_OVRSCAN2(ehu, &ovrscan2) != FH_SUCCESS ||
              fh_get_NAXIS(ehu, &naxis) != FH_SUCCESS ||
-             fh_get_NAXIS1(ehu, &naxis1) != FH_SUCCESS ||
-             fh_get_NAXIS2(ehu, &naxis2) != FH_SUCCESS) {
+             fh_get_NAXIS1(ehu, &n1) != FH_SUCCESS ||
+             fh_get_NAXIS2(ehu, &n2) != FH_SUCCESS) {
             fprintf(stderr, "error: Cannot get NAXIS* in extension %d row %d col %d\n", ext, row, col);
             exit(EXIT_FAILURE);
@@ -157 +187 @@
 /* Make some space now that we finally know the size */
          if(firstime) {
+            naxis1 = n1;
+            naxis2 = n2;
             firstime = 0;
             for(i=0; i<NCELL; i++) {
@@ -167 +199 @@
             Bbuf = (double *)calloc(naxis2, sizeof(double));
             Xbuf = (double *)calloc(NPOLY*naxis2, sizeof(double));
+         }
+         else {
+            if(naxis1 != n1 || naxis2 != n2) {
+               fprintf(stderr, "error: All cells must be same size (%dx%d != %dx%d)\n",
+                       n1, n2, naxis1, naxis2);
+               exit(EXIT_FAILURE);
+            }
          }
 
@@ -214 +253 @@
 #ifdef DEBUG
       for(col=0; col<NCELL; col++) {
-         printf("%1d %1d %3d %3d %6d %6d %6d %6d %8.3f\n",
+         printf("%1d %1d %3d %3d 0x%08x %6d %6d %6d %6d %8.3f\n",
                 col, row, otarow[col].cell, otarow[col].mefxt, 
+                (int)otarow[col].hu,
                 otarow[col].bias, otarow[col].sky, otarow[col].rms, 
                 otarow[col].nsat, otarow[col].ampl);
@@ -258 +298 @@
 /* Return to all the cells in the MEF and rewrite the data */
       for(col=0; col<NCELL; col++) {
+#ifdef DEBUG
+         printf("Now writing cell %d %d ext %d (hu 0x%08x)\n", 
+                col, row, otarow[col].mefxt, (int)otarow[col].hu);
+#endif
          if(otarow[col].mefxt < 0) continue;
+/* Skip if cell mask not set */
+         if(cellmask[col+NCELL*row] != 1) continue;
 /* This seeks back to the start of data from this cell again... */
          ehu = otarow[col].hu;
          fh_ehu(ehu, 0);
 
+#ifdef DEBUG
+         printf("Now writing cell %d %d ext %d (hu 0x%08x)\n", 
+                col, row, otarow[col].mefxt, (int)otarow[col].hu);
+#endif
          if (fh_write_padded_image(ehu, fh_file_desc(ehu), otarow[col].data,
                                    naxis1*naxis2*sizeof(short),
@@ -272 +322 @@
       }
 
-   }
+   }    /* Loop over rows */
 
    fh_destroy(ihu);
@@ -439 +489 @@
       A = Abuf[j] - (par[0] + j*par[1] + j*j*par[2]);
 #ifdef DEBUG
-      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f\n", j, Abuf[j], A);
+//      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f\n", j, Abuf[j], A);
 #endif
       for(k=0; k<NCELL; k++) {
@@ -504 +554 @@
       B = Bbuf[j] - (Bpar[0] + j*Bpar[1] + j*j*Bpar[2]);
 #ifdef DEBUG
-      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f %9.3f %9.3f\n", 
-                                      j, Abuf[j], A, Bbuf[j], B);
+//      if(otarow[3].mefxt == 4) printf("%4d %9.3f %9.3f %9.3f %9.3f\n", 
+//                                    j, Abuf[j], A, Bbuf[j], B);
 #endif
 /* Parabola passes through 0 at i=sx+nx, A at sx+nx/6, B at sx+nx/2 */
@@ -590 +640 @@
 
 #ifdef DEBUG
-      if(otarow[3].mefxt == 4) printf("%4d %4d %9.3f %9.3f %9.3f\n", 
-                                      j, i0, A, B, C);
+//      if(otarow[3].mefxt == 4) printf("%4d %4d %9.3f %9.3f %9.3f\n", 
+//                                    j, i0, A, B, C);
 #endif
 
@@ -609 +659 @@
    return(0);
 }
+
+static void syntax(const char *prog)
+{
+   printf("Syntax: %s mef_file [options]\n", prog);
+   printf("   where options include:\n\n");
+   printf(" mask=0101...   64 digits to work on cells 0:63.\n");
+   printf(" alg=N          Use fit algorithm 1=linear, 2=parabola (default)\n");   printf("                   3=bent linear (not yet fully implemented)\n");
+}