Index: /trunk/stac/src/Makefile =================================================================== --- /trunk/stac/src/Makefile (revision 3668) +++ /trunk/stac/src/Makefile (revision 3669) @@ -9,6 +9,5 @@ stacHelp.o -SHIFT = shift.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \ - stacCombine.o stacInvertMaps.o stacSize.o +SHIFT = shift.o stacRead.o stacTransform.o stacCheckMemory.o stacInvertMaps.o stacSize.o SHIFTSIZE = shiftSize.o stacWrite.o stacConfig.o stacRead.o stacCheckMemory.o stacSize.o Index: /trunk/stac/src/shift.c =================================================================== --- /trunk/stac/src/shift.c (revision 3668) +++ /trunk/stac/src/shift.c (revision 3669) @@ -1,13 +1,23 @@ #include +#include +#include +#include +#include #include "pslib.h" #include "stac.h" -#include -double getTime(void) -/* Gets the current time. Got this from Nick Kaiser's fetchpix.c */ +void help(const char *programName + ) { - struct timeval tv; - gettimeofday(&tv, NULL); - return(tv.tv_sec + 1.e-6 * tv.tv_usec); + fprintf (stderr, "shift: shift an image, given the transformation\n" + "Usage: %s [-h] [-v] [-s NX NY] OUT IN\n" + "where\n" + "\t-h Help (this info)\n" + "\t-v Increase verbosity level\n" + "\t-s NX NY Size of output image\n" + "\tOUT Output image\n" + "\tIN Input image, which has associated .map file.\n", + programName + ); } @@ -15,86 +25,118 @@ int main(int argc, char **argv) { - - double startTime = getTime(); - -#if 0 - psMemAllocateCallbackSet(stacMemPrint); - psMemAllocateCallbackSetID(185); - psMemFreeCallbackSet(stacMemPrint); - psMemFreeCallbackSetID(185); -#endif - // Set trace levels - psTraceSetLevel(".",0); - psTraceSetLevel("stac.checkMemory",10); - psTraceSetLevel("stac.config",10); - psTraceSetLevel("stac.read",10); - psTraceSetLevel("stac.invertMaps",10); - psTraceSetLevel("stac.errors",10); - psTraceSetLevel("stac.transform",10); - psTraceSetLevel("stac.combine",10); - psTraceSetLevel("stac.rejection",10); - psTraceSetLevel("stac.time",10); - psTraceSetLevel("stac.area",10); - psTraceSetLevel("stac.size",10); + psTraceSetLevel(".", 0); + psTraceSetLevel("stac.checkMemory", 10); + psTraceSetLevel("stac.read", 10); + psTraceSetLevel("stac.invertMaps", 10); + psTraceSetLevel("stac.transform", 10); + psTraceSetLevel("stac.size", 10); + psTraceSetLevel("shift", 10); // Set logging level psLogSetLevel(9); + // Parameters with default values + int outnx = 0, outny = 0; // Size of output image + int verbose = 0; // Verbosity level + // Parse command line - stacConfig *config = stacParseConfig(argc, argv); + 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 */ - // Read input files - psArray *inputs = stacReadImages(config); + /* Parse command-line arguments using getopt */ + while ((opt = getopt(argc, argv, "hvs:")) != -1) { + switch (opt) { + case 'h': + help(programName); + exit(EXIT_SUCCESS); + case 'v': + verbose++; + break; + 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. + help(programName); + exit(EXIT_FAILURE); + } + break; + default: + help(programName); + exit(EXIT_FAILURE); + } + } - // Read maps - psArray *maps = stacReadMaps(config); + /* Get the remaining, mandatory values */ + argc -= optind; + argv += optind; + if (argc < 2) { + help(programName); + exit(EXIT_FAILURE); + } + const char *inName = argv[0]; // Input filename + const char *outName = argv[1]; // Output filename - psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime); + psFits *imageFile = psFitsAlloc(inName); + psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read + psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0); + if (image == NULL) { + 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; + } + psFree(imageFile); + psFree(imageRegion); + + // Read map + 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 + images->data[0] = image; + maps->data[0] = map; // Get size - stacSize(config, inputs, maps); + if (outnx == 0 || outny == 0) { + stacSize(&outnx, &outny, NULL, NULL, images, maps); + } // Invert maps - psArray *inverseMaps = stacInvertMaps(maps, config); - - // Generate errors - psArray *errors = stacErrorImages(inputs, config); + psArray *inverseMaps = stacInvertMaps(maps, outnx, outny); // Transform inputs and errors psArray *transformed = NULL; - psArray *transformedErrors = NULL; - (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, config); + (void)stacTransform(&transformed, NULL, images, inverseMaps, NULL, NULL, NULL, NULL, NULL, outnx, outny); - psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime); - - // Combine with rejection - psArray *rejected = NULL; - psImage *combined = NULL; - (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, config); - - psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime); - - psFits *outFile = psFitsAlloc(config->output); - if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output); + // 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); } - psTrace("stac", 1, "Combined image written to %s\n", config->output); + psTrace("shift", 1, "Transformed image written to %s\n", outName); psFree(outFile); // Free everything I've used - stacConfigFree(config); - psFree(inputs); + psFree(images); psFree(maps); psFree(inverseMaps); - psFree(errors); - psFree(transformedErrors); psFree(transformed); - psFree(rejected); - psFree(combined); - psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime); +#ifdef TESTING + // Check memory for leaks, corruption + stacCheckMemory(); +#endif - // Check memory for leaks, corruption - stacCheckMemory(); } Index: /trunk/stac/src/stacSize.c =================================================================== --- /trunk/stac/src/stacSize.c (revision 3668) +++ /trunk/stac/src/stacSize.c (revision 3669) @@ -12,4 +12,5 @@ int nImages = images->n; // Number of input images assert(maps->n == nImages); + assert(outnx && outny); // Must be able to write to these psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates @@ -92,6 +93,10 @@ // Tweak the maps to account for the offset - *xMapDiff = floor(xMin); - *yMapDiff = floor(yMin); + if (xMapDiff != NULL) { + *xMapDiff = floor(xMin); + } + if (yMapDiff != NULL) { + *yMapDiff = floor(yMin); + } for (int i = 0; i < nImages; i++) { psPlaneTransform *map = maps->data[i]; // The map of interest Index: /trunk/stac/src/stacTransform.c =================================================================== --- /trunk/stac/src/stacTransform.c (revision 3668) +++ /trunk/stac/src/stacTransform.c (revision 3669) @@ -166,7 +166,11 @@ psImage *image = images->data[n]; // The input image psPlaneTransform *map = maps->data[n]; // The map - psImage *error = errors->data[n]; // The error image psImage *outImage = (*outputs)->data[n]; // The output image - psImage *outError = (*outErrors)->data[n]; // The output error image + psImage *error = NULL; // The error image + psImage *outError = NULL; // The output error image + if (errors) { + error = errors->data[n]; + outError = (*outErrors)->data[n]; + } float offset = 0.0; // Relative offset float scale = 1.0; // Relative scale @@ -177,14 +181,4 @@ scale = scales->data.F32[n]; } - -#if 0 - // No need for initialisation, since we iterate over the entire output image. - for (int y = 0; y < outny; y++) { - for (int x = 0; x < outnx; x++) { - outImage->data.F32[y][x] = 0.0; - outError->data.F32[y][x] = 0.0; - } - } -#endif // Mask @@ -208,9 +202,11 @@ detector->y, mask, 1, 0.0, PS_INTERPOLATE_BILINEAR); - // Error is actually the variance - outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, - detector->x, - detector->y, - mask, 1, 0.0); + if (error) { + // Error is actually the variance + outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, + detector->x, + detector->y, + mask, 1, 0.0); + } outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;