Index: trunk/stac/src/shift.c =================================================================== --- trunk/stac/src/shift.c (revision 3673) +++ trunk/stac/src/shift.c (revision 5743) @@ -11,13 +11,13 @@ { 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,9 +38,10 @@ // 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,13 +50,13 @@ /* 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,9 +64,15 @@ 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,39 +84,53 @@ 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,17 +140,29 @@ // 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);