IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3666


Ignore:
Timestamp:
Apr 5, 2005, 11:51:26 AM (21 years ago)
Author:
Paul Price
Message:

Removing stac-specific configuration out of functions so that I can use them in other programs. Program appears to work as it did before.

Location:
trunk/stac/src
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/stac.c

    r3660 r3666  
    22#include "pslib.h"
    33#include "stac.h"
     4#include "stacConfig.h"
    45#include <sys/time.h>
    56
    … …  
    4748
    4849    // Read input files
    49     psArray *inputs = stacReadImages(config);
     50    psArray *inputs = stacReadImages(config->inputs);
    5051
    5152    // Read maps
    52     psArray *maps = stacReadMaps(config);
     53    psArray *maps = stacReadMaps(config->inputs);
    5354
    5455    psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime);
    … …  
    5657    // Get size, if not input
    5758    if (config->outnx == 0 || config->outny == 0) {
    58         stacSize(config, inputs, maps);
     59        stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps);
    5960    }
    6061
    6162    // Invert maps
    62     psArray *inverseMaps = stacInvertMaps(maps, config);
     63    psArray *inverseMaps = stacInvertMaps(maps, config->outnx, config->outny);
    6364
    6465    // Generate errors
    65     psArray *errors = stacErrorImages(inputs, config);
     66    psArray *errors = stacErrorImages(inputs, config->gain, config->readnoise);
     67#ifdef TESTING
     68    // Write error images out to check
     69    for (int i = 0; i < inputs->n; i++) {
     70        char errName[MAXCHAR];          // Filename of error image
     71        sprintf(errName,"%s.err",config->inputs->data[i]);
     72        psFits *errorFile = psFitsAlloc(errName);
     73        if (!psFitsWriteImage(errorFile, NULL, error, 0, NULL)) {
     74            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
     75        }
     76        psTrace("stac", 1, "Error image written to %s\n", errName);
     77        psFree(errorFile);
     78    }
     79#endif
    6680
    6781    // Transform inputs and errors
    … …  
    6983    psArray *transformedErrors = NULL;
    7084    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL,
    71                         config);
    72 
     85                        config->outnx, config->outny);
    7386    psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime);
     87
     88#ifdef TESTING
     89    // Write transformed images out to check
     90    for (int i = 0; i < inputs->n; i++) {
     91        char shiftName[MAXCHAR];        // Filename of shift image
     92        char errName[MAXCHAR];          // Filename of error image
     93        sprintf(shiftName,"%s.shift1",config->inputs->data[n],numTransforms);
     94        sprintf(errName,"%s.shifterr1",config->inputs->data[n],numTransforms);
     95        psTrace("stac.transform.test", 6,
     96                "Output files have size: %dx%d\n",outImage->numCols,outImage->numRows);
     97        psFits *shiftFile = psFitsAlloc(shiftName);
     98        psFits *errFile = psFitsAlloc(errName);
     99        if (!psFitsWriteImage(shiftFile, NULL, outImage, 0, NULL)) {
     100            psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
     101        }
     102        psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
     103        if (!psFitsWriteImage(errFile, NULL, outError, 0, NULL)) {
     104            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
     105        }
     106        psTrace("stac", 1, "Shifted error image written to %s\n", errName);
     107        psFree(shiftFile);
     108        psFree(errFile);
     109    }
     110#endif
     111
     112
    74113
    75114    psVector *scales = NULL;            // Relative scales between images
    76115    psVector *offsets = NULL;           // Offsets between images
    77     (void)stacScales(&scales, &offsets, transformed, config);
     116    (void)stacScales(&scales, &offsets, transformed, config->starFile, config->starMapFile, config->xMapDiff,
     117                     config->yMapDiff, config->aper);
    78118    // Rescale the images
    79119    (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets);
     120    // Set the saturation and bad values
     121    psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits
     122    psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits
     123    for (int i = 0; i < transformed->n; i++) {
     124        saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
     125        bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
     126    }
    80127
    81128    // Combine with rejection
    82129    psArray *rejected = NULL;
    83130    psImage *combined = NULL;
    84     (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, config);
     131    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, saturated,
     132                      bad, config->reject);
    85133
    86134    psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime);
    87135
    88 #ifdef TESTING
     136
     137#ifdef TESTING
     138    // Write rejection images out to check
     139    for (int i = 0; i < nImages; i++) {
     140        char rejName[MAXCHAR];  // Filename of rejection image
     141        sprintf(rejName,"%s.shiftrej",config->inputs->data[i]);
     142       
     143        psFits *rejFile = psFitsAlloc(rejName);
     144        if (!psFitsWriteImage(rejFile, NULL, (*rejected)->data[i], 0, NULL)) {
     145            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
     146        }
     147        psTrace("stac", 1, "Rejection image written to %s\n", rejName);
     148        psFree(rejFile);
     149    }
     150
     151    // Write out pre-combined image
    89152    char preName[MAXCHAR];              // Filename of precombined image
    90153    sprintf(preName,"%s.pre",config->output);
    … …  
    116179
    117180    // Transform rejected pixels to source frame
    118     psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config);
     181    psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config->frac,
     182                                            config->grad, config->inputs);
    119183
    120184    // Get regions of interest in the output frame
    … …  
    134198    // Redo transformation with the masks and scales/offsets
    135199    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,
    136                         combineRegion, scales, offsets, config);
     200                        combineRegion, scales, offsets, config->outnx, config->outny);
    137201
    138202    // Combine the newly-transformed CR-free images, no rejection
    139203    psFree(rejected);
    140204    rejected = NULL;
    141     (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, config);
     205    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, saturated,
     206                      bad, config->reject);
    142207
    143208    // Write output image

  • trunk/stac/src/stac.h

    r3660 r3666  
    1313/************************************************************************************************************/
    1414
    15 // stacConfig.c
    16 
    17 // Configuration options
    18 typedef struct {
    19     int verbose;                        // Verbosity level
    20     float gain, readnoise;              // Gain and readnoise for detectors
    21     psArray *inputs;                    // Input file names
    22     const char *output;                 // Output file name
    23     bool saveShifts;                    // Save shifted images?
    24     const char *starFile;               // File with star coordinates
    25     const char *starMapFile;            // File with map for stars
    26     int outnx, outny;                   // Size of output image
    27     float xOffset, yOffset;             // Offset to get lower-left corner at 0,0
    28     psVector *saturated;                // Saturation level for each image
    29     psVector *bad;                      // Bad level for each image
    30     float reject;                       // Rejection level
    31     float frac;                         // Fraction of input pixel that must be masked before the pixel is
    32                                         // considered bad
    33     float grad;                         // Multiplier of the gradient
    34     int nReject;                        // Number of rejection iterations
    35     float xMapDiff, yMapDiff;           // Difference between pure map and output image coordinates
    36     float aper;                         // Aperture size (pixels)
    37 } stacConfig;
    38 
    39 // Allocator
    40 stacConfig *stacConfigAlloc(void);
    41 // Deallocator
    42 void stacConfigFree(stacConfig *config);
    43 
    44 // Help message
    45 void help(const char *name);
    46 
    47 // Parse the command line and return config
    48 stacConfig *stacParseConfig(int argc,   // Number of command-line arguments
    49                             char **argv // Command-line arguments
    50     );
    51 
    52 
    53 /************************************************************************************************************/
    54 
    5515// stacRead.c
    5616
    5717// Read the input files and return an array of images
    58 psArray *stacReadImages(stacConfig *config);
     18psArray *stacReadImages(psArray *filenames // The file names of the images
     19    );
    5920
    6021// Read a file of coordinates (x,y)
    … …  
    6526
    6627// Read the map files and return an array of transformations
    67 psArray *stacReadMaps(stacConfig *config);
     28psArray *stacReadMaps(psArray *filenames // The file names of the images whose maps are to be read
     29    );
    6830
    6931/************************************************************************************************************/
    … …  
    7335// Calculate the error images
    7436psArray *stacErrorImages(psArray *inputs, // Array of input images
    75                          stacConfig *config // Configuration details
     37                         float gain,    // Gain, in e/ADU
     38                         float rn       // Read noise, in e
    7639    );
    7740
    … …  
    9053                   const psVector *scales, // Relative scales
    9154                   const psVector *offsets, // Relative offsets
    92                    const stacConfig *config // Configuration
     55                   int outnx, int outny // Size of output images
    9356    );
    9457
    … …  
    13396                 int nReject,           // Number of rejection iterations
    13497                 psImage *region,       // Region to combine
    135                  stacConfig *config     // Configuration
     98                 psVector *saturated,   // Saturation limits for each image
     99                 psVector *bad,         // Bad pixel limits for each image
     100                 float reject           // Rejection (k-sigma)
    136101    );
    137102
    … …  
    142107// Invert an array of maps
    143108psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
    144                         const stacConfig *config // Configuration
    145     );
    146 
     109                        int outnx, int outny // Size of output image
     110    );
    147111
    148112/************************************************************************************************************/
    … …  
    161125                       psArray *maps,   // Maps from input to transformed image
    162126                       psArray *inverseMaps, // Maps from transformed to input image
    163                        stacConfig *config // Configuration
     127                       float frac,      // Fraction of pixel rejected before looking more carefully
     128                       float grad,      // Gradient limit for rejection
     129                       psArray *names   // Names of original images (only for writing out when TESTING)
    164130    );
    165131
    … …  
    186152
    187153// Calculate the size of the output image
    188 bool stacSize(stacConfig *config,       // Configuration, containing the output size
     154bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
     155              float *xMapDiff, float *yMapDiff, // Difference applied to maps
    189156              const psArray *images,    // Input images
    190157              psArray *maps             // Transformation maps
    191158    );
    192 
    193159
    194160/************************************************************************************************************/
    … …  
    206172                psVector **offsetsPtr,  // Offsets to return
    207173                const psArray *images,  // Images on which to measure the scales and offsets
    208                 const stacConfig *config // Configuration
     174                const char *starFile,   // File containing coordinates to photometer
     175                const char *starMapFile, // Map for coodinates to the common output frame
     176                float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
     177                float aper              // Aperture to use for photometry (radius)
    209178    );
    210179
    … …  
    224193
    225194// Write maps out
    226 bool stacWriteMaps(const psArray *maps, // Maps to write
    227                    const stacConfig *config // Configuration
     195bool stacWriteMap(const char *mapName,  // Filename to write to
     196                  psPlaneTransform *map // Map to write
    228197    );
    229198

  • trunk/stac/src/stacCombine.c

    r3610 r3666  
    6262                 int nReject,           // Number of rejection iterations
    6363                 psImage *region,       // Region to combine
    64                  stacConfig *config     // Configuration
     64                 psVector *saturated,   // Saturation limits for each image
     65                 psVector *bad,         // Bad pixel limits for each image
     66                 float reject           // Rejection (k-sigma)
    6567    )
    6668{
    … …  
    6870    int numRows = ((psImage*)images->data[0])->numRows; // Image size
    6971    int numCols = ((psImage*)images->data[0])->numCols; // Image size
    70     psVector *saturated = config->saturated;// Saturation limits
    71     psVector *bad = config->bad;        // Bad pixel limits
    72     float reject = config->reject;      // Rejection (k-sigma)
    7372   
    7473    // Check dimensions for consistency
    … …  
    188187
    189188#ifdef TESTING
    190     // Write rejection images out to check
    191     if (nReject > 0) {
    192         for (int i = 0; i < nImages; i++) {
    193             char rejName[MAXCHAR];      // Filename of rejection image
    194             sprintf(rejName,"%s.shiftrej",config->inputs->data[i]);
    195 
    196             psFits *rejFile = psFitsAlloc(rejName);
    197             if (!psFitsWriteImage(rejFile, NULL, (*rejected)->data[i], 0, NULL)) {
    198                 psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
    199             }
    200             psTrace("stac", 1, "Rejection image written to %s\n", rejName);
    201             psFree(rejFile);
    202         }
    203     }
    204189    // Write chi^2 image
    205190    iteration++;
    … …  
    210195        psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
    211196    }
    212     psTrace("stac", 1, "Chi^2 image written to %s\n", chiName);
     197    psTrace("stac.combine", 1, "Chi^2 image written to %s\n", chiName);
    213198    psFree(chiFile);
    214199    psFree(chi2);
    215200#endif
    216 
    217201
    218202    // Clean up

  • trunk/stac/src/stacConfig.c

    r3660 r3666  
    66#include "pslib.h"
    77#include "stac.h"
     8#include "stacConfig.h"
    89
    910stacConfig *stacConfigAlloc(void)
    … …  
    2122    config->starMapFile = NULL;
    2223    config->aper = 3.0;
    23     config->outnx = 512;
    24     config->outny = 512;
    25     config->saturated = NULL;           // Saturations; will fill this in later, when we know how many images
    26     config->bad = NULL;                 // Bad levels; will fill this in later, when we know how many images
     24    config->outnx = 0;
     25    config->outny = 0;
     26    config->saturated = 65535.0;        // Saturation level
     27    config->bad = 0.0;                  // Bad level
    2728    config->reject = 3.0;
    2829    config->frac = 0.5;
    … …  
    4041        psFree(config->inputs);
    4142    }
    42     if (config->saturated) {
    43         psFree(config->saturated);
    44     }
    45     if (config->bad) {
    46         psFree(config->bad);
    47     }
    4843    // Free everything
    4944    psFree(config);
    … …  
    5752    stacConfig *config = stacConfigAlloc(); // Configuration values
    5853    const char *programName = argv[0];  // Program name
    59 
    60     float saturated = 65535.0;          // Saturation level
    61     float bad = 0.0;                    // Bad level
    6254
    6355    /* Variables for getopt */
    … …  
    9688            break;
    9789          case 's':
    98             if (sscanf(optarg, "%f", &saturated) != 1) {
     90            if (sscanf(optarg, "%f", &config->saturated) != 1) {
    9991                help(programName);
    10092                exit(EXIT_FAILURE);
    … …  
    10294            break;
    10395          case 'b':
    104             if (sscanf(optarg, "%f", &bad) != 1) {
     96            if (sscanf(optarg, "%f", &config->bad) != 1) {
    10597                help(programName);
    10698                exit(EXIT_FAILURE);
    … …  
    170162    }
    171163
    172     // Create the saturation and bad vectors
    173     config->saturated = psVectorAlloc(argc-1, PS_TYPE_F32);
    174     config->bad = psVectorAlloc(argc-1, PS_TYPE_F32);
    175     for (int i = 0; i < argc - 1; i++) {
    176         ((psVector*)config->saturated)->data.F32[i] = saturated;
    177         ((psVector*)config->bad)->data.F32[i] = bad;
    178     }
    179 
    180164    // Debugging output
    181165    psTrace("stac.config", 8, "Parsed command line for configuration\n");

  • trunk/stac/src/stacErrorImages.c

    r3375 r3666  
    44
    55psArray *stacErrorImages(psArray *inputs, // Array of input images
    6                          stacConfig *config // Configuration details
     6                         float gain,    // Gain, in e/ADU
     7                         float rn       // Read noise, in e
    78    )
    89{
    9     float invGain = 1.0/config->gain;   // Inverse square root of gain
    10     float rn = config->readnoise/config->gain; // Read noise in ADU
     10    float invGain = 1.0/gain;           // Inverse square root of gain
     11    rn /= gain;                         // Read noise in ADU
    1112    psArray *errors = psArrayAlloc(inputs->n);
    1213
    … …  
    3233        // Put image onto the array
    3334        errors->data[i] = error;
    34  
    35 #ifdef TESTING
    36         // Write error image out to check
    37         char errName[MAXCHAR];          // Filename of error image
    38         sprintf(errName,"%s.err",config->inputs->data[i]);
    39         psFits *errorFile = psFitsAlloc(errName);
    40         if (!psFitsWriteImage(errorFile, NULL, error, 0, NULL)) {
    41             psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
    42         }
    43         psTrace("stac", 1, "Error image written to %s\n", errName);
    44         psFree(errorFile);
    45 #endif
    46 
    4735    }
    4836

  • trunk/stac/src/stacInvertMaps.c

    r3375 r3666  
    99
    1010psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
    11                         const stacConfig *config // Configuration
     11                        int outnx, int outny // Size of output image
    1212    )
    1313{
    … …  
    4848        // Create grid of points
    4949        for (int yint = 0; yint < NUM_GRID; yint++) {
    50             inCoord->y = (float)(yint * config->outny) / (float)(NUM_GRID - 1);
     50            inCoord->y = (float)(yint * outny) / (float)(NUM_GRID - 1);
    5151            for (int xint = 0; xint < NUM_GRID; xint++) {
    52                 inCoord->x = (float)(xint * config->outnx) / (float)(NUM_GRID - 1);
     52                inCoord->x = (float)(xint * outnx) / (float)(NUM_GRID - 1);
    5353
    5454                (void)psPlaneTransformApply(outCoord, oldMap, inCoord);

  • trunk/stac/src/stacRead.c

    r3610 r3666  
    66#define BUFFER 100                      // Size of buffer for incrementally reading coordinates
    77
    8 psArray *stacReadImages(stacConfig *config)
     8psArray *stacReadImages(psArray *filenames // The file names of the images
     9    )
    910{
    10     psArray *filenames = config->inputs;// The file names
    1111    int nFiles = filenames->n;          // The number of input files
    1212    psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
    … …  
    162162
    163163
    164 psArray *stacReadMaps(stacConfig *config)
     164psArray *stacReadMaps(psArray *filenames // The file names of the images whose maps are to be read
     165    )
    165166{
    166     psArray *filenames = config->inputs;// The file names
    167167    int nFiles = filenames->n;          // The number of input files
    168168    psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned

  • trunk/stac/src/stacRejection.c

    r3387 r3666  
    4848                       psArray *maps,   // Maps from input to transformed image
    4949                       psArray *inverseMaps, // Maps from transformed to input image
    50                        stacConfig *config // Configuration
     50                       float frac,      // Fraction of pixel rejected before looking more carefully
     51                       float grad,      // Gradient limit for rejection
     52                       psArray *names   // Names of original images (only for writing out when TESTING)
    5153    )
    5254{
    … …  
    5860    assert(inputs->n == maps->n);
    5961    assert(inputs->n == inverseMaps->n);
     62    assert(!names || names->n == inputs->n);
    6063
    6164    for (int i = 0; i < nImages; i++) {
    … …  
    9396        FILE *crs = NULL;               // File for outputting details of rejected pixels
    9497        char crfile[MAXCHAR];   // Filename
    95         sprintf(crfile,"%s.cr",config->inputs->data[i]);
     98        sprintf(crfile,"%s.cr",names->data[i]);
    9699        if ((crs = fopen(crfile, "w")) == NULL) {
    97100            fprintf(stderr, "Unable to open file for detailed output, %s\n");
    … …  
    119122#endif
    120123
    121                     if (maskVal > config->frac) {
     124                    if (maskVal > frac) {
    122125                        // Calculate mean gradient on other images
    123126                        float meanGrads = 0.0;
    … …  
    147150#endif
    148151
    149                         if (stacGradient(inputs->data[i], x, y) < config->grad * meanGrads) {
     152                        if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
    150153                            mask->data.U8[y][x] = 1;
    151154                            nBad++;
    … …  
    182185        char rejmapName[MAXCHAR];       // Filename of rejection image
    183186        char gradName[MAXCHAR];         // Filename of gradient image
    184         sprintf(maskName,"%s.mask",config->inputs->data[i]);
    185         sprintf(rejmapName,"%s.rejmap",config->inputs->data[i]);
    186         sprintf(gradName,"%s.grad",config->inputs->data[i]);
     187        sprintf(maskName, "%s.mask", names->data[i]);
     188        sprintf(rejmapName, "%s.rejmap", names->data[i]);
     189        sprintf(gradName, "%s.grad", names->data[i]);
    187190
    188191        psFits *maskFile = psFitsAlloc(maskName);

  • trunk/stac/src/stacScales.c

    r3613 r3666  
    5454                psVector **offsetsPtr,  // Offsets to return
    5555                const psArray *images,  // Images on which to measure the scales and offsets
    56                 const stacConfig *config // Configuration
     56                const char *starFile,   // File containing coordinates to photometer
     57                const char *starMapFile, // Map for coodinates to the common output frame
     58                float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
     59                float aper              // Aperture to use for photometry (radius)
    5760    )
    5861{
    … …  
    6063    assert(offsetsPtr);
    6164    assert(images);
    62     assert(config);
    6365    for (int i = 0; i < images->n; i++) {
    6466        psImage *image = images->data[i];
    … …  
    99101
    100102    // Now the scales
    101     if (config->starFile == NULL || config->starMapFile == NULL) {
     103    if (starFile == NULL || starMapFile == NULL) {
    102104        psLogMsg("stac.scales", PS_LOG_INFO,
    103105                 "No coordinates available to set scales --- assuming all are identical.\n");
    … …  
    108110    } else {
    109111        // Read star coordinates and map
    110         psArray *starCoords = stacReadCoords(config->starFile); // Array of star coordinates
    111         psPlaneTransform *starMap = stacReadMap(config->starMapFile); // Transformation for star coordinates
     112        psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates
     113        psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates
    112114
    113115        // Transform the stellar positions to match the transformed reference frame
    114116        psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
    115117        // Fix up difference between map and output frame
    116         starMap->x->coeff[0][0] -= config->xMapDiff;
    117         starMap->y->coeff[0][0] -= config->yMapDiff;
     118        starMap->x->coeff[0][0] -= xMapDiff;
     119        starMap->y->coeff[0][0] -= yMapDiff;
    118120        for (int i = 0; i < starCoords->n; i++) {
    119121            starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]);
    … …  
    139141                psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
    140142               
    141                 if (coords->x < config->aper || coords->y < config->aper ||
    142                     coords->x + config->aper > image->numCols - 1 ||
    143                     coords->y + config->aper > image->numRows) {
     143                if (coords->x < aper || coords->y < aper ||
     144                    coords->x + aper > image->numCols - 1 ||
     145                    coords->y + aper > image->numRows) {
    144146                    mask->data.U8[j] = 1;
    145147                } else {
    … …  
    147149                    float sum = 0.0;
    148150                    int numPix = 0;
    149                     float aper2 = SQUARE(config->aper);
    150                     for (int y = (int)floorf(coords->y - config->aper);
    151                          y <= (int)ceilf(coords->y + config->aper); y++) {
    152                         for (int x = (int)floorf(coords->x - config->aper);
    153                              x <= (int)ceilf(coords->x + config->aper); x++) {
     151                    float aper2 = SQUARE(aper);
     152                    for (int y = (int)floorf(coords->y - aper);
     153                         y <= (int)ceilf(coords->y + aper); y++) {
     154                        for (int x = (int)floorf(coords->x - aper);
     155                             x <= (int)ceilf(coords->x + aper); x++) {
    154156                            if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <=
    155157                                aper2) {
    … …  
    198200    }
    199201
    200     // Change the saturation and bad values
    201     psVector *saturated = config->saturated; // Saturation limits
    202     psVector *bad = config->bad;        // Bad limits
    203     for (int i = 0; i < saturated->n; i++) {
    204         saturated->data.F32[i] = (saturated->data.F32[i] - offsets->data.F32[i]) / scales->data.F32[i];
    205         bad->data.F32[i] = (bad->data.F32[i] - offsets->data.F32[i]) / scales->data.F32[i];
    206     }
    207 
    208202    return true;
    209203}

  • trunk/stac/src/stacSize.c

    r3610 r3666  
    44#include "stac.h"
    55
    6 bool stacSize(stacConfig *config,       // Configuration, containing the output size
     6bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
     7              float *xMapDiff, float *yMapDiff, // Difference applied to maps
    78              const psArray *images,    // Input images
    89              psArray *maps             // Transformation maps
    … …  
    9192
    9293    // Tweak the maps to account for the offset
    93     config->xMapDiff = floor(xMin);
    94     config->yMapDiff = floor(yMin);
     94    *xMapDiff = floor(xMin);
     95    *yMapDiff = floor(yMin);
    9596    for (int i = 0; i < nImages; i++) {
    9697        psPlaneTransform *map = maps->data[i]; // The map of interest
    … …  
    99100    }
    100101
    101     config->outnx = (int)(xMax + 0.5) - (int)xMin;
    102     config->outny = (int)(yMax + 0.5) - (int)yMin;
     102    *outnx = (int)(xMax + 0.5) - (int)xMin;
     103    *outny = (int)(yMax + 0.5) - (int)yMin;
    103104
    104     psTrace("stac.size", 1, "Output size is to be %dx%d\n", config->outnx, config->outny);
     105    psTrace("stac.size", 1, "Output size is to be %dx%d\n", *outnx, *outny);
    105106
    106107    psFree(inCoord);

  • trunk/stac/src/stacTransform.c

    r3610 r3666  
    77#define MIN(x,y) (((x) > (y)) ? (y) : (x))
    88#define MAX(x,y) (((x) > (y)) ? (x) : (y))
    9 
    10 #define MAXCHAR 80
    11 
    12 static int numTransforms = 0;           // Number of transformations performed
    13 
    149
    1510// Hacked the original ps_ImagePixelInterpolateBILINEAR_F32 to add variances
    … …  
    115110                   const psVector *scales, // Relative scales
    116111                   const psVector *offsets, // Relative offsets
    117                    const stacConfig *config // Configuration
     112                   int outnx, int outny // Size of output images
    118113    )
    119114{
    120115    int nImages = images->n;            // Number of images
    121     int nx = config->outnx, ny = config->outny; // Size of output images
    122     numTransforms++;
    123116
    124117    // Check input sizes
    … …  
    134127    if (*outputs == NULL) {
    135128        *outputs = psArrayAlloc(nImages);
    136         psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", nx, ny);
     129        psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
    137130        for (int i = 0; i < nImages; i++) {
    138             (*outputs)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
     131            (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
    139132        }
    140133    }
    … …  
    144137    if (errors && (*outErrors == NULL)) {
    145138        *outErrors = psArrayAlloc(errors->n);
    146         psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", nx, ny);
     139        psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
    147140        for (int i = 0; i < nImages; i++) {
    148             (*outErrors)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
     141            (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
    149142        }
    150143    }
    … …  
    187180#if 0
    188181        // No need for initialisation, since we iterate over the entire output image.
    189         for (int y = 0; y < ny; y++) {
    190             for (int x = 0; x < nx; x++) {
     182        for (int y = 0; y < outny; y++) {
     183            for (int x = 0; x < outnx; x++) {
    191184                outImage->data.F32[y][x] = 0.0;
    192185                outError->data.F32[y][x] = 0.0;
    … …  
    202195
    203196        // Iterate over the output image pixels
    204         for (int y = 0; y < ny; y++) {
    205             for (int x = 0; x < nx; x++) {
     197        for (int y = 0; y < outny; y++) {
     198            for (int x = 0; x < outnx; x++) {
    206199                // Only transform those pixels requested
    207200                if (!region || (region && region->data.U8[y][x])) {
    … …  
    229222        } // Iterating over output pixels
    230223
    231 #ifdef TESTING
    232         // Write error image out to check
    233         char shiftName[MAXCHAR];        // Filename of shift image
    234         char errName[MAXCHAR];          // Filename of error image
    235         sprintf(shiftName,"%s.shift.%d",config->inputs->data[n],numTransforms);
    236         sprintf(errName,"%s.shifterr.%d",config->inputs->data[n],numTransforms);
    237         psTrace("stac.transform.test", 6,
    238                 "Output files have size: %dx%d\n",outImage->numCols,outImage->numRows);
    239 
    240         psFits *shiftFile = psFitsAlloc(shiftName);
    241         psFits *errFile = psFitsAlloc(errName);
    242         if (!psFitsWriteImage(shiftFile, NULL, outImage, 0, NULL)) {
    243             psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
    244         }
    245         psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
    246         if (!psFitsWriteImage(errFile, NULL, outError, 0, NULL)) {
    247             psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
    248         }
    249         psTrace("stac", 1, "Shifted error image written to %s\n", errName);
    250         psFree(shiftFile);
    251         psFree(errFile);
    252 #endif
    253 
    254224    } // Iterating over images
    255225
Note: See TracChangeset for help on using the changeset viewer.