Index: trunk/stac/src/combine.c =================================================================== --- trunk/stac/src/combine.c (revision 3680) +++ trunk/stac/src/combine.c (revision 3680) @@ -0,0 +1,70 @@ +#include +#include "pslib.h" +#include "stac.h" +#include "combineConfig.h" + +int main(int argc, char **argv) +{ + // Set trace levels + psTraceSetLevel(".", 0); + psTraceSetLevel("stac.checkMemory", 10); + psTraceSetLevel("stac.read", 10); + psTraceSetLevel("stac.scales", 7); + psTraceSetLevel("stac.rescale", 10); + psTraceSetLevel("stac.combine", 10); + psTraceSetLevel("combine", 10); + + // Set logging level + psLogSetLevel(9); + + // Parse command line + combineConfig *config = combineParseConfig(argc, argv); // Configuration + + psArray *images = stacReadImages(config->inNames); // The images + psArray *errors = stacErrorImages(images, config->gain, config->readnoise); // Error images + + // Calculate scales between images + psVector *scales = NULL; // Relative scales between images + psVector *offsets = NULL; // Offsets between images + (void)stacScales(&scales, &offsets, images, config->starFile, config->starMap, 0.0, 0.0, config->aper); + + // Set the saturation and bad values + psVector *saturated = psVectorAlloc(images->n, PS_TYPE_F32); // Saturation limits + psVector *bad = psVectorAlloc(images->n, PS_TYPE_F32); // Bad limits + for (int i = 0; i < images->n; i++) { + saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i]; + bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i]; + } + + // Rescale the images + (void)stacRescale(images, errors, NULL, scales, offsets); + + // Combine the images + psImage *combined = NULL; // Combined image + psArray *rejected = NULL; // Array of rejection masks + stacCombine(&combined, &rejected, images, errors, config->nReject, NULL, saturated, bad, config->reject); + + psFits *outFile = psFitsAlloc(config->outName); + if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", config->outName); + } + psTrace("combine", 1, "Combined image written to %s\n", config->outName); + psFree(outFile); + + // Clean up + psFree(combined); + psFree(rejected); + psFree(images); + psFree(errors); + psFree(saturated); + psFree(bad); + psFree(scales); + psFree(offsets); + combineConfigFree(config); + +#ifdef TESTING + // Check memory for leaks, corruption + stacCheckMemory(); +#endif + +} Index: trunk/stac/src/combineConfig.c =================================================================== --- trunk/stac/src/combineConfig.c (revision 3680) +++ trunk/stac/src/combineConfig.c (revision 3680) @@ -0,0 +1,161 @@ +#include +#include +#include +#include +#include +#include "pslib.h" +#include "combineConfig.h" + +void help(const char *programName + ) +{ + fprintf (stderr, "shift: shift an image, given the transformation\n" + "Usage: %s [-h] [-v] [-g GAIN] [-r READNOISE] [-s SAT] [-b BAD] [-p FILE MAP] [-a APER] [-k SIGMAREJ] [-n NREJECT] OUT IN1 IN2...\n" + "where\n" + "\t-h Help (this info)\n" + "\t-v Increase verbosity level\n" + "\t-g GAIN Gain in e/ADU (1.0)\n" + "\t-r READNOISE Read noise in e (0.0)\n" + "\t-s SAT Saturation point (65535)\n" + "\t-b BAD Bad level (0)\n" + "\t-p FILE MAP Specify file containing star coordinates, with map\n" + "\t-a APER Aperture radius for photometry (3.0)\n" + "\t-k SIGMAREJ k-sigma rejection threshold (3.0)\n" + "\t-n NREJECT Number of rejection iterations (1)\n" + "\tOUT Output image\n" + "\tIN1 IN2... Input images (identical size)\n", + programName + ); +} + + +combineConfig *combineConfigAlloc(void) +{ + combineConfig *config = psAlloc(sizeof(combineConfig)); // Configuration + + // Parameters with default values + config->verbose = 0; // Verbosity level + config->gain = 1.0; // Gain (e/ADU) + config->readnoise = 0.0; // Read noise (e) + config->reject = 4.0; // Rejection threshold (sigma) + config->nReject = 1; // Number of rejection iterations + config->saturated = 65535.0; // Saturation level + config->bad = 0.0; // Bad level + config->outName = NULL; // Output name + config->inNames = NULL; // Input names; + config->starFile = NULL; // Filename of file containing stars + config->starMap = NULL; // Map for stars + config->aper = 3.0; // Aperture for photometry + + return config; +} + +void combineConfigFree(combineConfig *config) +{ + psFree(config->inNames); + psFree(config); +} + + +combineConfig *combineParseConfig(int argc, char *argv[]) +{ + combineConfig *config = combineConfigAlloc(); // Configuration + + const char *programName = argv[0]; // Program name + + /* Variables for getopt */ + int opt; /* Option, from getopt */ + extern char *optarg; /* Argument accompanying switch */ + extern int optind; /* getopt variables */ + + /* Parse command-line arguments using getopt */ + while ((opt = getopt(argc, argv, "hvg:r:s:b:p:a:k:n:")) != -1) { + switch (opt) { + case 'h': + help(programName); + exit(EXIT_SUCCESS); + case 'v': + config->verbose++; + break; + case 'r': + if (sscanf(optarg, "%f", &config->readnoise) != 1) { + printf("Unable to read readnoise.\n"); + help(programName); + exit(EXIT_FAILURE); + } + break; + case 'g': + if (sscanf(optarg, "%f", &config->gain) != 1) { + printf("Unable to read gain.\n"); + help(programName); + exit(EXIT_FAILURE); + } + break; + case 's': + if (sscanf(optarg, "%f", &config->saturated) != 1) { + printf("Unable to read saturation limit.\n"); + help(programName); + exit(EXIT_FAILURE); + } + break; + case 'b': + if (sscanf(optarg, "%f", &config->bad) != 1) { + printf("Unable to read bad limit.\n"); + help(programName); + exit(EXIT_FAILURE); + } + break; + case 'p': + if (argc < optind+1) { + printf("Unable to read photometric files.\n"); + help(programName); + exit(EXIT_FAILURE); + } + config->starFile = argv[optind-1]; + config->starMap = argv[optind++]; + // Note: incrementing optind, so I can read more than one parameter. + break; + case 'a': + if (sscanf(optarg, "%f", &config->aper) != 1) { + printf("Unable to read aperture.\n"); + help(programName); + exit(EXIT_FAILURE); + } + break; + case 'k': + if (sscanf(optarg, "%f", &config->reject) != 1) { + printf("Unable to read rejection limit.\n"); + help(programName); + exit(EXIT_FAILURE); + } + break; + case 'n': + if (sscanf(optarg, "%d", &config->nReject) != 1) { + printf("Unable to read number of rejection iterations.\n"); + help(programName); + exit(EXIT_FAILURE); + } + break; + default: + printf("Bad option: %c\n", opt); + help(programName); + exit(EXIT_FAILURE); + } + } + + /* Get the remaining, mandatory values */ + argc -= optind; + argv += optind; + if (argc < 2) { + help(programName); + exit(EXIT_FAILURE); + } + config->outName = argv[0]; // Output filename + config->inNames = psArrayAlloc(argc-1); // Input filenames + for (int i = 1; i < argc; i++) { + config->inNames->data[i-1] = psAlloc(strlen(argv[i])); + strncpy(config->inNames->data[i-1], argv[i], strlen(argv[i])); + } + + return config; +} Index: trunk/stac/src/combineConfig.h =================================================================== --- trunk/stac/src/combineConfig.h (revision 3680) +++ trunk/stac/src/combineConfig.h (revision 3680) @@ -0,0 +1,28 @@ +#ifndef COMBINE_CONFIG_H +#define COMBINE_CONFIG_H + +typedef struct { + int verbose; // Verbosity level + float gain; // Gain (e/ADU) + float readnoise; // Read noise (e) + float reject; // Rejection threshold (sigma) + int nReject; // Number of rejection iterations + float saturated; // Saturation level + float bad; // Bad level + char *outName; // Output name + psArray *inNames; // Input names; + char *starFile; // Filename of file containing stars + char *starMap; // Map for stars + float aper; // Aperture for photometry +} combineConfig; + + +combineConfig *combineConfigAlloc(void);// Allocator +void combineConfigFree(combineConfig *config); // Deallocator + +void help(const char *programName); // Print help +combineConfig *combineParseConfig(int argc, char *argv[]); // Parse command line + + + +#endif Index: trunk/stac/src/stacRead.c =================================================================== --- trunk/stac/src/stacRead.c (revision 3673) +++ trunk/stac/src/stacRead.c (revision 3680) @@ -32,4 +32,8 @@ psFree(image); } else { + int numNaN = psImageClipNaN(image, 0.0); + if (numNaN) { + psTrace("stac.read.images", 5, "Clipped %d NaN pixels.\n", numNaN); + } images->data[i] = image; } Index: trunk/stac/src/stacWrite.c =================================================================== --- trunk/stac/src/stacWrite.c (revision 3680) +++ trunk/stac/src/stacWrite.c (revision 3680) @@ -0,0 +1,77 @@ +#include +#include +#include "pslib.h" +#include "stac.h" + +bool stacWriteMap(const char *mapName, // Filename to write to + psPlaneTransform *map // Map to write + ) +{ + assert(mapName); + + FILE *mapFile = fopen(mapName, "w"); + if (!mapFile) { + fprintf(stderr, "Unable to open map file: %s\n", mapName); + return false; + } + + psDPolynomial2D *xMap = map->x; // x transform + psDPolynomial2D *yMap = map->y; // y transform + + // A crucial limitation of the current system --- the order of each polynomial must be the same + assert(xMap->nX == xMap->nY && yMap->nX == yMap->nY && xMap->nX == yMap->nX); + int order = xMap->nX - 1; // The polynomial order + fprintf(mapFile, "%d\n", order); + + // x coefficients + for (int k = 0; k < order + 1; k++) { + for (int j = 0; j < k + 1; j++) { + int i = k - j; + if (xMap->mask[i][j]) { + fprintf(mapFile, "0.0 "); + } else { + fprintf(mapFile, "%g ", xMap->coeff[i][j]); + } + } + } + fprintf(mapFile, "\n"); + + // y coefficients + for (int k = 0; k < order + 1; k++) { + for (int j = 0; j < k + 1; j++) { + int i = k - j; + if (yMap->mask[i][j]) { + fprintf(mapFile, "0.0 "); + } else { + fprintf(mapFile, "%g ", yMap->coeff[i][j]); + } + } + } + fprintf(mapFile, "\n"); + + fclose(mapFile); + + return true; +} + + +bool stacWriteMaps(const psArray *names, // Filenames of the input images (will add ".map") + const psArray *maps // Maps to write + ) +{ + assert(names); + assert(maps); + assert(names->n == maps->n); + + for (int i = 0; i < names->n; i++) { + char mapName[MAXCHAR]; // Filename of error image + sprintf(mapName, "%s.map", names->data[i]); + if (!stacWriteMap(mapName, maps->data[i])) { + return false; + } + } + + return true; +} + +