Changeset 5743 for trunk/stac/src/shift.c

Timestamp:

Dec 7, 2005, 2:57:14 PM (20 years ago)

Author:

Paul Price

Message:

Importing from PAP cvs tree again. This will now be the working tree.

File:

• trunk/stac/src/shift.c (modified) (6 diffs)

trunk/stac/src/shift.c

@@ -11 +11 @@
 {
     fprintf (stderr, "shift: shift an image, given the transformation\n"
-             "Usage: %s [-h] [-v] [-s NX NY] IN OUT\n"
-             "where\n"
-             "\t-h           Help (this info)\n"
-             "\t-v           Increase verbosity level\n"
-             "\t-s NX NY     Size of output image\n"
-             "\tIN           Input image, which has associated .map file\n"
-             "\tOUT          Output image\n",
-             programName
-        );
+             "Usage: %s [-h] [-v] [-s NX NY] IN OUT\n"
+             "where\n"
+             "\t-h           Help (this info)\n"
+             "\t-v           Increase verbosity level\n"
+             "\t-s NX NY     Size of output image\n"
+             "\tIN           Input image, which has associated .map file\n"
+             "\tOUT          Output image\n",
+             programName
+        );
 }
 
@@ -38 +38 @@
 
     // Parameters with default values
-    int outnx = 0, outny = 0;           // Size of output image
-    int verbose = 0;                    // Verbosity level
+    int outnx = 0, outny = 0;           // Size of output image
+    int verbose = 0;                    // Verbosity level
+    float bad = 0.0;                    // Bad level
 
     // Parse command line
-    const char *programName = argv[0];  // Program name
+    const char *programName = argv[0];  // Program name
     /* Variables for getopt */
     int opt;   /* Option, from getopt */
@@ -49 +50 @@
 
     /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvs:")) != -1) {
+    while ((opt = getopt(argc, argv, "hvb:s:")) != -1) {
         switch (opt) {
-          case 'h':
-            help(programName);
-            exit(EXIT_SUCCESS);
-          case 'v':
+          case 'h':
+            help(programName);
+            exit(EXIT_SUCCESS);
+          case 'v':
             verbose++;
             break;
-          case 's':
+          case 's':
             if ((sscanf(argv[optind-1], "%d", &outnx) != 1) || (sscanf(argv[optind++], "%d", &outny) != 1)) {
                 // Note: incrementing optind, so I can read more than one parameter.
@@ -63 +64 @@
                 exit(EXIT_FAILURE);
             }
-            break;
-          default:
-            help(programName);
-            exit(EXIT_FAILURE);
-        }
+            break;
+          case 'b':
+            if (sscanf(optarg, "%f", &bad) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          default:
+            help(programName);
+            exit(EXIT_FAILURE);
+        }
     }
 
@@ -77 +84 @@
         exit(EXIT_FAILURE);
     }
-    const char *inName = argv[0];       // Input filename
-    const char *outName = argv[1];      // Output filename
+    const char *inName = argv[0];       // Input filename
+    const char *outName = argv[1];      // Output filename
 
-    psFits *imageFile = psFitsAlloc(inName);
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
-    psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
+    psFits *imageFile = psFitsOpen(inName, "r");
+    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
+    psMetadata *header = psFitsReadHeader(NULL, imageFile); // FITS header
+    psImage *image = psFitsReadImage(NULL, imageFile, imageRegion, 0);
     if (image == NULL) {
-        psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", inName);
-        exit(EXIT_FAILURE);
+        psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", inName);
+        exit(EXIT_FAILURE);
     }
     psTrace("shift", 4, "Image %s is %dx%d\n", inName, image->numCols, image->numRows);
     // Convert to 32-bit floating point, in necessary
     if (image->type.type != PS_TYPE_F32) {
-        psTrace("shift", 5, "Converting %s to floating point in memory....", inName);
-        psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
-        psFree(image);
-        image = temp;
+        psTrace("shift", 5, "Converting %s to floating point in memory....", inName);
+        psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
+        psFree(image);
+        image = temp;
     }
-    psFree(imageFile);
-    psFree(imageRegion);
+    psFitsClose(imageFile);
+
+    // Generate masks
+    psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
+    for (int y = 0; y < image->numRows; y++) {
+        for (int x = 0; x < image->numCols; x++) {
+            if (image->data.F32[y][x] <= bad) {
+                mask->data.U8[y][x] = 1;
+            } else {
+                mask->data.U8[y][x] = 0;
+            }
+        }
+    }
 
     // Read map
-    char mapName[MAXCHAR];              // Name of map file
+    char mapName[MAXCHAR];              // Name of map file
     sprintf(mapName, "%s.map", inName);
     psPlaneTransform *map = stacReadMap(mapName);
 
     // Functions work on array of images, so:
-    psArray *images = psArrayAlloc(1); // An array of one
-    psArray *maps = psArrayAlloc(1);// An array of one
+    psArray *images = psArrayAlloc(1);  // An array of one
+    psArray *maps = psArrayAlloc(1);    // An array of one
+    psArray *masks = psArrayAlloc(1);   // An array of one
     images->data[0] = image;
     maps->data[0] = map;
+    masks->data[0] = mask;
 
     // Get size
     if (outnx == 0 || outny == 0) {
-        stacSize(&outnx, &outny, NULL, NULL, images, maps);
+        stacSize(&outnx, &outny, NULL, NULL, images, maps);
     }
 
@@ -119 +140 @@
     // Transform inputs and errors
     psArray *transformed = NULL;
-    (void)stacTransform(&transformed, NULL, images, inverseMaps, NULL, NULL, NULL, NULL, NULL, outnx, outny);
+    (void)stacTransform(&transformed, NULL, images, inverseMaps, NULL, masks, NULL, NULL, NULL, outnx, outny);
+
+    // Update FITS header appropriately
+    psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of pixels in x",
+                     (int)((psImage*)transformed->data[0])->numCols);
+    psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of pixels in y",
+                     (int)((psImage*)transformed->data[0])->numRows);
+    psMetadataAddS32(header, PS_LIST_TAIL, "BITPIX", PS_META_REPLACE, "Bits per pixel", -32);
 
     // Write out transformed image
-    psFits *outFile = psFitsAlloc(outName);
-    if (!psFitsWriteImage(outFile, NULL, transformed->data[0], 0, NULL)) {
-        psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
+    psFits *outFile = psFitsOpen(outName, "w");
+    int numPix = psImageClipNaN(transformed->data[0], 0.0);
+    if (numPix > 0) {
+        psTrace("stac", 3, "Clipping %d NaN pixels to zero.\n", numPix);
+    }
+    if (!psFitsWriteImage(outFile, header, transformed->data[0], 0)) {
+        psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
     }
     psTrace("shift", 1, "Transformed image written to %s\n", outName);
-    psFree(outFile);
+    psFitsClose(outFile);
 
     // Free everything I've used
     psFree(images);
     psFree(maps);
+    psFree(masks);
     psFree(inverseMaps);
     psFree(transformed);