Index: /trunk/stac/src/stac.c =================================================================== --- /trunk/stac/src/stac.c (revision 3665) +++ /trunk/stac/src/stac.c (revision 3666) @@ -2,4 +2,5 @@ #include "pslib.h" #include "stac.h" +#include "stacConfig.h" #include @@ -47,8 +48,8 @@ // Read input files - psArray *inputs = stacReadImages(config); + psArray *inputs = stacReadImages(config->inputs); // Read maps - psArray *maps = stacReadMaps(config); + psArray *maps = stacReadMaps(config->inputs); psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime); @@ -56,12 +57,25 @@ // Get size, if not input if (config->outnx == 0 || config->outny == 0) { - stacSize(config, inputs, maps); + stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps); } // Invert maps - psArray *inverseMaps = stacInvertMaps(maps, config); + psArray *inverseMaps = stacInvertMaps(maps, config->outnx, config->outny); // Generate errors - psArray *errors = stacErrorImages(inputs, config); + psArray *errors = stacErrorImages(inputs, config->gain, config->readnoise); +#ifdef TESTING + // Write error images out to check + for (int i = 0; i < inputs->n; i++) { + char errName[MAXCHAR]; // Filename of error image + sprintf(errName,"%s.err",config->inputs->data[i]); + psFits *errorFile = psFitsAlloc(errName); + if (!psFitsWriteImage(errorFile, NULL, error, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); + } + psTrace("stac", 1, "Error image written to %s\n", errName); + psFree(errorFile); + } +#endif // Transform inputs and errors @@ -69,22 +83,71 @@ psArray *transformedErrors = NULL; (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL, - config); - + config->outnx, config->outny); psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime); + +#ifdef TESTING + // Write transformed images out to check + for (int i = 0; i < inputs->n; i++) { + char shiftName[MAXCHAR]; // Filename of shift image + char errName[MAXCHAR]; // Filename of error image + sprintf(shiftName,"%s.shift1",config->inputs->data[n],numTransforms); + sprintf(errName,"%s.shifterr1",config->inputs->data[n],numTransforms); + psTrace("stac.transform.test", 6, + "Output files have size: %dx%d\n",outImage->numCols,outImage->numRows); + psFits *shiftFile = psFitsAlloc(shiftName); + psFits *errFile = psFitsAlloc(errName); + if (!psFitsWriteImage(shiftFile, NULL, outImage, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); + } + psTrace("stac", 1, "Shifted image written to %s\n", shiftName); + if (!psFitsWriteImage(errFile, NULL, outError, 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); + } + psTrace("stac", 1, "Shifted error image written to %s\n", errName); + psFree(shiftFile); + psFree(errFile); + } +#endif + + psVector *scales = NULL; // Relative scales between images psVector *offsets = NULL; // Offsets between images - (void)stacScales(&scales, &offsets, transformed, config); + (void)stacScales(&scales, &offsets, transformed, config->starFile, config->starMapFile, config->xMapDiff, + config->yMapDiff, config->aper); // Rescale the images (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets); + // Set the saturation and bad values + psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits + psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits + for (int i = 0; i < transformed->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]; + } // Combine with rejection psArray *rejected = NULL; psImage *combined = NULL; - (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, config); + (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, saturated, + bad, config->reject); psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime); -#ifdef TESTING + +#ifdef TESTING + // Write rejection images out to check + for (int i = 0; i < nImages; i++) { + char rejName[MAXCHAR]; // Filename of rejection image + sprintf(rejName,"%s.shiftrej",config->inputs->data[i]); + + psFits *rejFile = psFitsAlloc(rejName); + if (!psFitsWriteImage(rejFile, NULL, (*rejected)->data[i], 0, NULL)) { + psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName); + } + psTrace("stac", 1, "Rejection image written to %s\n", rejName); + psFree(rejFile); + } + + // Write out pre-combined image char preName[MAXCHAR]; // Filename of precombined image sprintf(preName,"%s.pre",config->output); @@ -116,5 +179,6 @@ // Transform rejected pixels to source frame - psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config); + psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config->frac, + config->grad, config->inputs); // Get regions of interest in the output frame @@ -134,10 +198,11 @@ // Redo transformation with the masks and scales/offsets (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource, - combineRegion, scales, offsets, config); + combineRegion, scales, offsets, config->outnx, config->outny); // Combine the newly-transformed CR-free images, no rejection psFree(rejected); rejected = NULL; - (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, config); + (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, saturated, + bad, config->reject); // Write output image Index: /trunk/stac/src/stac.h =================================================================== --- /trunk/stac/src/stac.h (revision 3665) +++ /trunk/stac/src/stac.h (revision 3666) @@ -13,48 +13,9 @@ /************************************************************************************************************/ -// stacConfig.c - -// Configuration options -typedef struct { - int verbose; // Verbosity level - float gain, readnoise; // Gain and readnoise for detectors - psArray *inputs; // Input file names - const char *output; // Output file name - bool saveShifts; // Save shifted images? - const char *starFile; // File with star coordinates - const char *starMapFile; // File with map for stars - int outnx, outny; // Size of output image - float xOffset, yOffset; // Offset to get lower-left corner at 0,0 - psVector *saturated; // Saturation level for each image - psVector *bad; // Bad level for each image - float reject; // Rejection level - float frac; // Fraction of input pixel that must be masked before the pixel is - // considered bad - float grad; // Multiplier of the gradient - int nReject; // Number of rejection iterations - float xMapDiff, yMapDiff; // Difference between pure map and output image coordinates - float aper; // Aperture size (pixels) -} stacConfig; - -// Allocator -stacConfig *stacConfigAlloc(void); -// Deallocator -void stacConfigFree(stacConfig *config); - -// Help message -void help(const char *name); - -// Parse the command line and return config -stacConfig *stacParseConfig(int argc, // Number of command-line arguments - char **argv // Command-line arguments - ); - - -/************************************************************************************************************/ - // stacRead.c // Read the input files and return an array of images -psArray *stacReadImages(stacConfig *config); +psArray *stacReadImages(psArray *filenames // The file names of the images + ); // Read a file of coordinates (x,y) @@ -65,5 +26,6 @@ // Read the map files and return an array of transformations -psArray *stacReadMaps(stacConfig *config); +psArray *stacReadMaps(psArray *filenames // The file names of the images whose maps are to be read + ); /************************************************************************************************************/ @@ -73,5 +35,6 @@ // Calculate the error images psArray *stacErrorImages(psArray *inputs, // Array of input images - stacConfig *config // Configuration details + float gain, // Gain, in e/ADU + float rn // Read noise, in e ); @@ -90,5 +53,5 @@ const psVector *scales, // Relative scales const psVector *offsets, // Relative offsets - const stacConfig *config // Configuration + int outnx, int outny // Size of output images ); @@ -133,5 +96,7 @@ int nReject, // Number of rejection iterations psImage *region, // Region to combine - stacConfig *config // Configuration + psVector *saturated, // Saturation limits for each image + psVector *bad, // Bad pixel limits for each image + float reject // Rejection (k-sigma) ); @@ -142,7 +107,6 @@ // Invert an array of maps psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert - const stacConfig *config // Configuration - ); - + int outnx, int outny // Size of output image + ); /************************************************************************************************************/ @@ -161,5 +125,7 @@ psArray *maps, // Maps from input to transformed image psArray *inverseMaps, // Maps from transformed to input image - stacConfig *config // Configuration + float frac, // Fraction of pixel rejected before looking more carefully + float grad, // Gradient limit for rejection + psArray *names // Names of original images (only for writing out when TESTING) ); @@ -186,9 +152,9 @@ // Calculate the size of the output image -bool stacSize(stacConfig *config, // Configuration, containing the output size +bool stacSize(int *outnx, int *outny, // Size of output image (returned) + float *xMapDiff, float *yMapDiff, // Difference applied to maps const psArray *images, // Input images psArray *maps // Transformation maps ); - /************************************************************************************************************/ @@ -206,5 +172,8 @@ psVector **offsetsPtr, // Offsets to return const psArray *images, // Images on which to measure the scales and offsets - const stacConfig *config // Configuration + const char *starFile, // File containing coordinates to photometer + const char *starMapFile, // Map for coodinates to the common output frame + float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference) + float aper // Aperture to use for photometry (radius) ); @@ -224,6 +193,6 @@ // Write maps out -bool stacWriteMaps(const psArray *maps, // Maps to write - const stacConfig *config // Configuration +bool stacWriteMap(const char *mapName, // Filename to write to + psPlaneTransform *map // Map to write ); Index: /trunk/stac/src/stacCombine.c =================================================================== --- /trunk/stac/src/stacCombine.c (revision 3665) +++ /trunk/stac/src/stacCombine.c (revision 3666) @@ -62,5 +62,7 @@ int nReject, // Number of rejection iterations psImage *region, // Region to combine - stacConfig *config // Configuration + psVector *saturated, // Saturation limits for each image + psVector *bad, // Bad pixel limits for each image + float reject // Rejection (k-sigma) ) { @@ -68,7 +70,4 @@ int numRows = ((psImage*)images->data[0])->numRows; // Image size int numCols = ((psImage*)images->data[0])->numCols; // Image size - psVector *saturated = config->saturated;// Saturation limits - psVector *bad = config->bad; // Bad pixel limits - float reject = config->reject; // Rejection (k-sigma) // Check dimensions for consistency @@ -188,18 +187,4 @@ #ifdef TESTING - // Write rejection images out to check - if (nReject > 0) { - for (int i = 0; i < nImages; i++) { - char rejName[MAXCHAR]; // Filename of rejection image - sprintf(rejName,"%s.shiftrej",config->inputs->data[i]); - - psFits *rejFile = psFitsAlloc(rejName); - if (!psFitsWriteImage(rejFile, NULL, (*rejected)->data[i], 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName); - } - psTrace("stac", 1, "Rejection image written to %s\n", rejName); - psFree(rejFile); - } - } // Write chi^2 image iteration++; @@ -210,9 +195,8 @@ psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName); } - psTrace("stac", 1, "Chi^2 image written to %s\n", chiName); + psTrace("stac.combine", 1, "Chi^2 image written to %s\n", chiName); psFree(chiFile); psFree(chi2); #endif - // Clean up Index: /trunk/stac/src/stacConfig.c =================================================================== --- /trunk/stac/src/stacConfig.c (revision 3665) +++ /trunk/stac/src/stacConfig.c (revision 3666) @@ -6,4 +6,5 @@ #include "pslib.h" #include "stac.h" +#include "stacConfig.h" stacConfig *stacConfigAlloc(void) @@ -21,8 +22,8 @@ config->starMapFile = NULL; config->aper = 3.0; - config->outnx = 512; - config->outny = 512; - config->saturated = NULL; // Saturations; will fill this in later, when we know how many images - config->bad = NULL; // Bad levels; will fill this in later, when we know how many images + config->outnx = 0; + config->outny = 0; + config->saturated = 65535.0; // Saturation level + config->bad = 0.0; // Bad level config->reject = 3.0; config->frac = 0.5; @@ -40,10 +41,4 @@ psFree(config->inputs); } - if (config->saturated) { - psFree(config->saturated); - } - if (config->bad) { - psFree(config->bad); - } // Free everything psFree(config); @@ -57,7 +52,4 @@ stacConfig *config = stacConfigAlloc(); // Configuration values const char *programName = argv[0]; // Program name - - float saturated = 65535.0; // Saturation level - float bad = 0.0; // Bad level /* Variables for getopt */ @@ -96,5 +88,5 @@ break; case 's': - if (sscanf(optarg, "%f", &saturated) != 1) { + if (sscanf(optarg, "%f", &config->saturated) != 1) { help(programName); exit(EXIT_FAILURE); @@ -102,5 +94,5 @@ break; case 'b': - if (sscanf(optarg, "%f", &bad) != 1) { + if (sscanf(optarg, "%f", &config->bad) != 1) { help(programName); exit(EXIT_FAILURE); @@ -170,12 +162,4 @@ } - // Create the saturation and bad vectors - config->saturated = psVectorAlloc(argc-1, PS_TYPE_F32); - config->bad = psVectorAlloc(argc-1, PS_TYPE_F32); - for (int i = 0; i < argc - 1; i++) { - ((psVector*)config->saturated)->data.F32[i] = saturated; - ((psVector*)config->bad)->data.F32[i] = bad; - } - // Debugging output psTrace("stac.config", 8, "Parsed command line for configuration\n"); Index: /trunk/stac/src/stacConfig.h =================================================================== --- /trunk/stac/src/stacConfig.h (revision 3666) +++ /trunk/stac/src/stacConfig.h (revision 3666) @@ -0,0 +1,43 @@ +#ifndef STAC_CONFIG_H +#define STAC_CONFIG_H + +#include "pslib.h" + +// stacConfig.c + +// Configuration options +typedef struct { + int verbose; // Verbosity level + float gain, readnoise; // Gain and readnoise for detectors + psArray *inputs; // Input file names + const char *output; // Output file name + bool saveShifts; // Save shifted images? + const char *starFile; // File with star coordinates + const char *starMapFile; // File with map for stars + int outnx, outny; // Size of output image + float xOffset, yOffset; // Offset to get lower-left corner at 0,0 + float saturated; // Saturation level + float bad; // Bad level + float reject; // Rejection level + float frac; // Fraction of input pixel that must be masked before the pixel is + // considered bad + float grad; // Multiplier of the gradient + int nReject; // Number of rejection iterations + float xMapDiff, yMapDiff; // Difference between pure map and output image coordinates + float aper; // Aperture size (pixels) +} stacConfig; + +// Allocator +stacConfig *stacConfigAlloc(void); +// Deallocator +void stacConfigFree(stacConfig *config); + +// Help message +void help(const char *name); + +// Parse the command line and return config +stacConfig *stacParseConfig(int argc, // Number of command-line arguments + char **argv // Command-line arguments + ); + +#endif Index: /trunk/stac/src/stacErrorImages.c =================================================================== --- /trunk/stac/src/stacErrorImages.c (revision 3665) +++ /trunk/stac/src/stacErrorImages.c (revision 3666) @@ -4,9 +4,10 @@ psArray *stacErrorImages(psArray *inputs, // Array of input images - stacConfig *config // Configuration details + float gain, // Gain, in e/ADU + float rn // Read noise, in e ) { - float invGain = 1.0/config->gain; // Inverse square root of gain - float rn = config->readnoise/config->gain; // Read noise in ADU + float invGain = 1.0/gain; // Inverse square root of gain + rn /= gain; // Read noise in ADU psArray *errors = psArrayAlloc(inputs->n); @@ -32,17 +33,4 @@ // Put image onto the array errors->data[i] = error; - -#ifdef TESTING - // Write error image out to check - char errName[MAXCHAR]; // Filename of error image - sprintf(errName,"%s.err",config->inputs->data[i]); - psFits *errorFile = psFitsAlloc(errName); - if (!psFitsWriteImage(errorFile, NULL, error, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); - } - psTrace("stac", 1, "Error image written to %s\n", errName); - psFree(errorFile); -#endif - } Index: /trunk/stac/src/stacInvertMaps.c =================================================================== --- /trunk/stac/src/stacInvertMaps.c (revision 3665) +++ /trunk/stac/src/stacInvertMaps.c (revision 3666) @@ -9,5 +9,5 @@ psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert - const stacConfig *config // Configuration + int outnx, int outny // Size of output image ) { @@ -48,7 +48,7 @@ // Create grid of points for (int yint = 0; yint < NUM_GRID; yint++) { - inCoord->y = (float)(yint * config->outny) / (float)(NUM_GRID - 1); + inCoord->y = (float)(yint * outny) / (float)(NUM_GRID - 1); for (int xint = 0; xint < NUM_GRID; xint++) { - inCoord->x = (float)(xint * config->outnx) / (float)(NUM_GRID - 1); + inCoord->x = (float)(xint * outnx) / (float)(NUM_GRID - 1); (void)psPlaneTransformApply(outCoord, oldMap, inCoord); Index: /trunk/stac/src/stacRead.c =================================================================== --- /trunk/stac/src/stacRead.c (revision 3665) +++ /trunk/stac/src/stacRead.c (revision 3666) @@ -6,7 +6,7 @@ #define BUFFER 100 // Size of buffer for incrementally reading coordinates -psArray *stacReadImages(stacConfig *config) +psArray *stacReadImages(psArray *filenames // The file names of the images + ) { - psArray *filenames = config->inputs;// The file names int nFiles = filenames->n; // The number of input files psArray *images = psArrayAlloc(nFiles); // The input files, to be returned @@ -162,7 +162,7 @@ -psArray *stacReadMaps(stacConfig *config) +psArray *stacReadMaps(psArray *filenames // The file names of the images whose maps are to be read + ) { - psArray *filenames = config->inputs;// The file names int nFiles = filenames->n; // The number of input files psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned Index: /trunk/stac/src/stacRejection.c =================================================================== --- /trunk/stac/src/stacRejection.c (revision 3665) +++ /trunk/stac/src/stacRejection.c (revision 3666) @@ -48,5 +48,7 @@ psArray *maps, // Maps from input to transformed image psArray *inverseMaps, // Maps from transformed to input image - stacConfig *config // Configuration + float frac, // Fraction of pixel rejected before looking more carefully + float grad, // Gradient limit for rejection + psArray *names // Names of original images (only for writing out when TESTING) ) { @@ -58,4 +60,5 @@ assert(inputs->n == maps->n); assert(inputs->n == inverseMaps->n); + assert(!names || names->n == inputs->n); for (int i = 0; i < nImages; i++) { @@ -93,5 +96,5 @@ FILE *crs = NULL; // File for outputting details of rejected pixels char crfile[MAXCHAR]; // Filename - sprintf(crfile,"%s.cr",config->inputs->data[i]); + sprintf(crfile,"%s.cr",names->data[i]); if ((crs = fopen(crfile, "w")) == NULL) { fprintf(stderr, "Unable to open file for detailed output, %s\n"); @@ -119,5 +122,5 @@ #endif - if (maskVal > config->frac) { + if (maskVal > frac) { // Calculate mean gradient on other images float meanGrads = 0.0; @@ -147,5 +150,5 @@ #endif - if (stacGradient(inputs->data[i], x, y) < config->grad * meanGrads) { + if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) { mask->data.U8[y][x] = 1; nBad++; @@ -182,7 +185,7 @@ char rejmapName[MAXCHAR]; // Filename of rejection image char gradName[MAXCHAR]; // Filename of gradient image - sprintf(maskName,"%s.mask",config->inputs->data[i]); - sprintf(rejmapName,"%s.rejmap",config->inputs->data[i]); - sprintf(gradName,"%s.grad",config->inputs->data[i]); + sprintf(maskName, "%s.mask", names->data[i]); + sprintf(rejmapName, "%s.rejmap", names->data[i]); + sprintf(gradName, "%s.grad", names->data[i]); psFits *maskFile = psFitsAlloc(maskName); Index: /trunk/stac/src/stacScales.c =================================================================== --- /trunk/stac/src/stacScales.c (revision 3665) +++ /trunk/stac/src/stacScales.c (revision 3666) @@ -54,5 +54,8 @@ psVector **offsetsPtr, // Offsets to return const psArray *images, // Images on which to measure the scales and offsets - const stacConfig *config // Configuration + const char *starFile, // File containing coordinates to photometer + const char *starMapFile, // Map for coodinates to the common output frame + float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference) + float aper // Aperture to use for photometry (radius) ) { @@ -60,5 +63,4 @@ assert(offsetsPtr); assert(images); - assert(config); for (int i = 0; i < images->n; i++) { psImage *image = images->data[i]; @@ -99,5 +101,5 @@ // Now the scales - if (config->starFile == NULL || config->starMapFile == NULL) { + if (starFile == NULL || starMapFile == NULL) { psLogMsg("stac.scales", PS_LOG_INFO, "No coordinates available to set scales --- assuming all are identical.\n"); @@ -108,12 +110,12 @@ } else { // Read star coordinates and map - psArray *starCoords = stacReadCoords(config->starFile); // Array of star coordinates - psPlaneTransform *starMap = stacReadMap(config->starMapFile); // Transformation for star coordinates + psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates + psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates // Transform the stellar positions to match the transformed reference frame psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions // Fix up difference between map and output frame - starMap->x->coeff[0][0] -= config->xMapDiff; - starMap->y->coeff[0][0] -= config->yMapDiff; + starMap->x->coeff[0][0] -= xMapDiff; + starMap->y->coeff[0][0] -= yMapDiff; for (int i = 0; i < starCoords->n; i++) { starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]); @@ -139,7 +141,7 @@ psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star - if (coords->x < config->aper || coords->y < config->aper || - coords->x + config->aper > image->numCols - 1 || - coords->y + config->aper > image->numRows) { + if (coords->x < aper || coords->y < aper || + coords->x + aper > image->numCols - 1 || + coords->y + aper > image->numRows) { mask->data.U8[j] = 1; } else { @@ -147,9 +149,9 @@ float sum = 0.0; int numPix = 0; - float aper2 = SQUARE(config->aper); - for (int y = (int)floorf(coords->y - config->aper); - y <= (int)ceilf(coords->y + config->aper); y++) { - for (int x = (int)floorf(coords->x - config->aper); - x <= (int)ceilf(coords->x + config->aper); x++) { + float aper2 = SQUARE(aper); + for (int y = (int)floorf(coords->y - aper); + y <= (int)ceilf(coords->y + aper); y++) { + for (int x = (int)floorf(coords->x - aper); + x <= (int)ceilf(coords->x + aper); x++) { if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <= aper2) { @@ -198,12 +200,4 @@ } - // Change the saturation and bad values - psVector *saturated = config->saturated; // Saturation limits - psVector *bad = config->bad; // Bad limits - for (int i = 0; i < saturated->n; i++) { - saturated->data.F32[i] = (saturated->data.F32[i] - offsets->data.F32[i]) / scales->data.F32[i]; - bad->data.F32[i] = (bad->data.F32[i] - offsets->data.F32[i]) / scales->data.F32[i]; - } - return true; } Index: /trunk/stac/src/stacSize.c =================================================================== --- /trunk/stac/src/stacSize.c (revision 3665) +++ /trunk/stac/src/stacSize.c (revision 3666) @@ -4,5 +4,6 @@ #include "stac.h" -bool stacSize(stacConfig *config, // Configuration, containing the output size +bool stacSize(int *outnx, int *outny, // Size of output image (returned) + float *xMapDiff, float *yMapDiff, // Difference applied to maps const psArray *images, // Input images psArray *maps // Transformation maps @@ -91,6 +92,6 @@ // Tweak the maps to account for the offset - config->xMapDiff = floor(xMin); - config->yMapDiff = floor(yMin); + *xMapDiff = floor(xMin); + *yMapDiff = floor(yMin); for (int i = 0; i < nImages; i++) { psPlaneTransform *map = maps->data[i]; // The map of interest @@ -99,8 +100,8 @@ } - config->outnx = (int)(xMax + 0.5) - (int)xMin; - config->outny = (int)(yMax + 0.5) - (int)yMin; + *outnx = (int)(xMax + 0.5) - (int)xMin; + *outny = (int)(yMax + 0.5) - (int)yMin; - psTrace("stac.size", 1, "Output size is to be %dx%d\n", config->outnx, config->outny); + psTrace("stac.size", 1, "Output size is to be %dx%d\n", *outnx, *outny); psFree(inCoord); Index: /trunk/stac/src/stacTransform.c =================================================================== --- /trunk/stac/src/stacTransform.c (revision 3665) +++ /trunk/stac/src/stacTransform.c (revision 3666) @@ -7,9 +7,4 @@ #define MIN(x,y) (((x) > (y)) ? (y) : (x)) #define MAX(x,y) (((x) > (y)) ? (x) : (y)) - -#define MAXCHAR 80 - -static int numTransforms = 0; // Number of transformations performed - // Hacked the original ps_ImagePixelInterpolateBILINEAR_F32 to add variances @@ -115,10 +110,8 @@ const psVector *scales, // Relative scales const psVector *offsets, // Relative offsets - const stacConfig *config // Configuration + int outnx, int outny // Size of output images ) { int nImages = images->n; // Number of images - int nx = config->outnx, ny = config->outny; // Size of output images - numTransforms++; // Check input sizes @@ -134,7 +127,7 @@ if (*outputs == NULL) { *outputs = psArrayAlloc(nImages); - psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", nx, ny); + psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny); for (int i = 0; i < nImages; i++) { - (*outputs)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32); + (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32); } } @@ -144,7 +137,7 @@ if (errors && (*outErrors == NULL)) { *outErrors = psArrayAlloc(errors->n); - psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", nx, ny); + psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny); for (int i = 0; i < nImages; i++) { - (*outErrors)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32); + (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32); } } @@ -187,6 +180,6 @@ #if 0 // No need for initialisation, since we iterate over the entire output image. - for (int y = 0; y < ny; y++) { - for (int x = 0; x < nx; x++) { + 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; @@ -202,6 +195,6 @@ // Iterate over the output image pixels - for (int y = 0; y < ny; y++) { - for (int x = 0; x < nx; x++) { + for (int y = 0; y < outny; y++) { + for (int x = 0; x < outnx; x++) { // Only transform those pixels requested if (!region || (region && region->data.U8[y][x])) { @@ -229,27 +222,4 @@ } // Iterating over output pixels -#ifdef TESTING - // Write error image out to check - char shiftName[MAXCHAR]; // Filename of shift image - char errName[MAXCHAR]; // Filename of error image - sprintf(shiftName,"%s.shift.%d",config->inputs->data[n],numTransforms); - sprintf(errName,"%s.shifterr.%d",config->inputs->data[n],numTransforms); - psTrace("stac.transform.test", 6, - "Output files have size: %dx%d\n",outImage->numCols,outImage->numRows); - - psFits *shiftFile = psFitsAlloc(shiftName); - psFits *errFile = psFitsAlloc(errName); - if (!psFitsWriteImage(shiftFile, NULL, outImage, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); - } - psTrace("stac", 1, "Shifted image written to %s\n", shiftName); - if (!psFitsWriteImage(errFile, NULL, outError, 0, NULL)) { - psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); - } - psTrace("stac", 1, "Shifted error image written to %s\n", errName); - psFree(shiftFile); - psFree(errFile); -#endif - } // Iterating over images