IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6887


Ignore:
Timestamp:
Apr 18, 2006, 12:20:45 PM (20 years ago)
Author:
Paul Price
Message:

Updating to use psLib rel11. Main problem was that vector and array lengths ('n' element) are no longer set to equal the number of allocated values ('nalloc' element) when allocated.

Location:
trunk/stac
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/configure.ac

    r6453 r6887  
    2626
    2727stac_CFLAGS="-Wall -Werror -std=c99 -DPS_NO_TRACE"
     28#stac_CFLAGS="-Wall -Werror -std=c99"
    2829AC_SUBST([stac_CFLAGS])
    2930

  • trunk/stac/src/combine.c

    r6771 r6887  
    2525    // Make fake errors
    2626    psArray *errors = psArrayAlloc(images->n);
     27    errors->n = images->n;
    2728    for (int i = 0; i < images->n; i++) {
    2829        psImage *image = images->data[i];
    … …  
    4445    psVector *saturated = psVectorAlloc(images->n, PS_TYPE_F32); // Saturation limits
    4546    psVector *bad = psVectorAlloc(images->n, PS_TYPE_F32); // Bad limits
     47    saturated->n = images->n;
     48    bad->n = images->n;
    4649    for (int i = 0; i < images->n; i++) {
    4750        saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];

  • trunk/stac/src/combineConfig.c

    r3680 r6887  
    1111{
    1212    fprintf (stderr, "shift: shift an image, given the transformation\n"
    13              "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"
    14              "where\n"
    15              "\t-h           Help (this info)\n"
    16              "\t-v           Increase verbosity level\n"
    17              "\t-g GAIN      Gain in e/ADU (1.0)\n"
    18              "\t-r READNOISE Read noise in e (0.0)\n"
    19              "\t-s SAT       Saturation point (65535)\n"
    20              "\t-b BAD       Bad level (0)\n"
    21              "\t-p FILE MAP  Specify file containing star coordinates, with map\n"
    22              "\t-a APER      Aperture radius for photometry (3.0)\n"
    23              "\t-k SIGMAREJ  k-sigma rejection threshold (3.0)\n"
    24              "\t-n NREJECT   Number of rejection iterations (1)\n"
    25              "\tOUT          Output image\n"
    26              "\tIN1 IN2...   Input images (identical size)\n",
    27              programName
    28         );
     13             "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"
     14             "where\n"
     15             "\t-h           Help (this info)\n"
     16             "\t-v           Increase verbosity level\n"
     17             "\t-g GAIN      Gain in e/ADU (1.0)\n"
     18             "\t-r READNOISE Read noise in e (0.0)\n"
     19             "\t-s SAT       Saturation point (65535)\n"
     20             "\t-b BAD       Bad level (0)\n"
     21             "\t-p FILE MAP  Specify file containing star coordinates, with map\n"
     22             "\t-a APER      Aperture radius for photometry (3.0)\n"
     23             "\t-k SIGMAREJ  k-sigma rejection threshold (3.0)\n"
     24             "\t-n NREJECT   Number of rejection iterations (1)\n"
     25             "\tOUT          Output image\n"
     26             "\tIN1 IN2...   Input images (identical size)\n",
     27             programName
     28        );
    2929}
    3030
    … …  
    3535
    3636    // Parameters with default values
    37     config->verbose = 0;                // Verbosity level
    38     config->gain = 1.0;                 // Gain (e/ADU)
    39     config->readnoise = 0.0;            // Read noise (e)
    40     config->reject = 4.0;               // Rejection threshold (sigma)
    41     config->nReject = 1;                // Number of rejection iterations
    42     config->saturated = 65535.0;        // Saturation level
    43     config->bad = 0.0;                  // Bad level
    44     config->outName = NULL;             // Output name
    45     config->inNames = NULL;             // Input names;
    46     config->starFile = NULL;            // Filename of file containing stars
    47     config->starMap = NULL;             // Map for stars
    48     config->aper = 3.0;                 // Aperture for photometry
     37    config->verbose = 0;                // Verbosity level
     38    config->gain = 1.0;                 // Gain (e/ADU)
     39    config->readnoise = 0.0;            // Read noise (e)
     40    config->reject = 4.0;               // Rejection threshold (sigma)
     41    config->nReject = 1;                // Number of rejection iterations
     42    config->saturated = 65535.0;        // Saturation level
     43    config->bad = 0.0;                  // Bad level
     44    config->outName = NULL;             // Output name
     45    config->inNames = NULL;             // Input names;
     46    config->starFile = NULL;            // Filename of file containing stars
     47    config->starMap = NULL;             // Map for stars
     48    config->aper = 3.0;                 // Aperture for photometry
    4949
    5050    return config;
    … …  
    6262    combineConfig *config = combineConfigAlloc(); // Configuration
    6363
    64     const char *programName = argv[0];  // Program name
     64    const char *programName = argv[0];  // Program name
    6565
    6666    /* Variables for getopt */
    … …  
    7272    while ((opt = getopt(argc, argv, "hvg:r:s:b:p:a:k:n:")) != -1) {
    7373        switch (opt) {
    74           case 'h':
    75             help(programName);
    76             exit(EXIT_SUCCESS);
    77           case 'v':
     74          case 'h':
     75            help(programName);
     76            exit(EXIT_SUCCESS);
     77          case 'v':
    7878            config->verbose++;
    7979            break;
    80           case 'r':
     80          case 'r':
    8181            if (sscanf(optarg, "%f", &config->readnoise) != 1) {
    82                 printf("Unable to read readnoise.\n");
    83                 help(programName);
    84                 exit(EXIT_FAILURE);
    85             }
    86             break;
    87           case 'g':
    88             if (sscanf(optarg, "%f", &config->gain) != 1) {
    89                 printf("Unable to read gain.\n");
    90                 help(programName);
    91                 exit(EXIT_FAILURE);
    92             }
    93             break;
    94           case 's':
    95             if (sscanf(optarg, "%f", &config->saturated) != 1) {
    96                 printf("Unable to read saturation limit.\n");
    97                 help(programName);
    98                 exit(EXIT_FAILURE);
    99             }
    100             break;
    101           case 'b':
    102             if (sscanf(optarg, "%f", &config->bad) != 1) {
    103                 printf("Unable to read bad limit.\n");
    104                 help(programName);
    105                 exit(EXIT_FAILURE);
    106             }
    107             break;
    108           case 'p':
    109             if (argc < optind+1) {
    110                 printf("Unable to read photometric files.\n");
    111                 help(programName);
    112                 exit(EXIT_FAILURE);
    113             }
    114             config->starFile = argv[optind-1];
    115             config->starMap = argv[optind++];
    116             // Note: incrementing optind, so I can read more than one parameter.
    117             break;
    118           case 'a':
    119             if (sscanf(optarg, "%f", &config->aper) != 1) {
    120                 printf("Unable to read aperture.\n");
     82                printf("Unable to read readnoise.\n");
    12183                help(programName);
    12284                exit(EXIT_FAILURE);
    12385            }
    124             break;
    125           case 'k':
    126             if (sscanf(optarg, "%f", &config->reject) != 1) {
    127                 printf("Unable to read rejection limit.\n");
     86            break;
     87          case 'g':
     88            if (sscanf(optarg, "%f", &config->gain) != 1) {
     89                printf("Unable to read gain.\n");
    12890                help(programName);
    12991                exit(EXIT_FAILURE);
    13092            }
    131             break;
    132           case 'n':
    133             if (sscanf(optarg, "%d", &config->nReject) != 1) {
    134                 printf("Unable to read number of rejection iterations.\n");
    135                 help(programName);
     93            break;
     94          case 's':
     95            if (sscanf(optarg, "%f", &config->saturated) != 1) {
     96                printf("Unable to read saturation limit.\n");
     97                help(programName);
    13698                exit(EXIT_FAILURE);
    13799            }
    138             break;
    139           default:
    140             printf("Bad option: %c\n", opt);
    141             help(programName);
    142             exit(EXIT_FAILURE);
    143         }
     100            break;
     101          case 'b':
     102            if (sscanf(optarg, "%f", &config->bad) != 1) {
     103                printf("Unable to read bad limit.\n");
     104                help(programName);
     105                exit(EXIT_FAILURE);
     106            }
     107            break;
     108          case 'p':
     109            if (argc < optind+1) {
     110                printf("Unable to read photometric files.\n");
     111                help(programName);
     112                exit(EXIT_FAILURE);
     113            }
     114            config->starFile = argv[optind-1];
     115            config->starMap = argv[optind++];
     116            // Note: incrementing optind, so I can read more than one parameter.
     117            break;
     118          case 'a':
     119            if (sscanf(optarg, "%f", &config->aper) != 1) {
     120                printf("Unable to read aperture.\n");
     121                help(programName);
     122                exit(EXIT_FAILURE);
     123            }
     124            break;
     125          case 'k':
     126            if (sscanf(optarg, "%f", &config->reject) != 1) {
     127                printf("Unable to read rejection limit.\n");
     128                help(programName);
     129                exit(EXIT_FAILURE);
     130            }
     131            break;
     132          case 'n':
     133            if (sscanf(optarg, "%d", &config->nReject) != 1) {
     134                printf("Unable to read number of rejection iterations.\n");
     135                help(programName);
     136                exit(EXIT_FAILURE);
     137            }
     138            break;
     139          default:
     140            printf("Bad option: %c\n", opt);
     141            help(programName);
     142            exit(EXIT_FAILURE);
     143        }
    144144    }
    145145
    … …  
    151151        exit(EXIT_FAILURE);
    152152    }
    153     config->outName = argv[0];          // Output filename
     153    config->outName = argv[0];          // Output filename
    154154    config->inNames = psArrayAlloc(argc-1); // Input filenames
     155    config->inNames->n = argc-1;
    155156    for (int i = 1; i < argc; i++) {
    156         config->inNames->data[i-1] = psAlloc(strlen(argv[i]));
    157         strncpy(config->inNames->data[i-1], argv[i], strlen(argv[i]));
     157        config->inNames->data[i-1] = psAlloc(strlen(argv[i]));
     158        strncpy(config->inNames->data[i-1], argv[i], strlen(argv[i]));
    158159    }
    159160

  • trunk/stac/src/shift.c

    r6771 r6887  
    126126    psArray *maps = psArrayAlloc(1);    // An array of one
    127127    psArray *masks = psArrayAlloc(1);   // An array of one
     128    images->n = maps->n = masks->n = 1;
    128129    images->data[0] = image;
    129130    maps->data[0] = map;
    … …  
    132133    // Get size
    133134    if (outnx == 0 || outny == 0) {
    134         stacSize(&outnx, &outny, NULL, NULL, images, maps);
     135        psVector *xSize = psVectorAlloc(1, PS_TYPE_S32); // A vector of one
     136        psVector *ySize = psVectorAlloc(1, PS_TYPE_S32); // A vector of one
     137        xSize->n = 1;
     138        ySize->n = 1;
     139        xSize->data.S32[0] = image->numCols;
     140        ySize->data.S32[0] = image->numRows;
     141        stacSize(&outnx, &outny, NULL, NULL, xSize, ySize, maps);
     142        psFree(xSize);
     143        psFree(ySize);
    135144    }
    136145

  • trunk/stac/src/shiftSize.c

    r5743 r6887  
    1717{
    1818    fprintf (stderr, "shiftSize: Calculate size for output combined image\n"
    19              "Usage: %s [-h] [-v] IN1 IN2...\n"
    20              "where\n"
    21              "\t-h           Help (this info)\n"
    22              "\t-v           Increase verbosity level\n"
    23              "\tIN1, IN2...  Input images, which have associated .map files.\n",
    24              name
    25         );
     19             "Usage: %s [-h] [-v] IN1 IN2...\n"
     20             "where\n"
     21             "\t-h           Help (this info)\n"
     22             "\t-v           Increase verbosity level\n"
     23             "\tIN1, IN2...  Input images, which have associated .map files.\n",
     24             name
     25        );
    2626}
    2727
    … …  
    4141
    4242
    43     int verbose = 0;                    // Verbosity level
     43    int verbose = 0;                    // Verbosity level
    4444
    4545    // Parse command line
    46     const char *programName = argv[0];  // Program name
     46    const char *programName = argv[0];  // Program name
    4747    /* Variables for getopt */
    4848    int opt;   /* Option, from getopt */
    … …  
    5353    while ((opt = getopt(argc, argv, "hv")) != -1) {
    5454        switch (opt) {
    55           case 'h':
    56             help(programName);
    57             exit(EXIT_SUCCESS);
    58           case 'v':
     55          case 'h':
     56            help(programName);
     57            exit(EXIT_SUCCESS);
     58          case 'v':
    5959            verbose++;
    6060            break;
    61           default:
    62             help(programName);
    63         }
     61          default:
     62            help(programName);
     63        }
    6464    }
    6565
    … …  
    7373
    7474    psArray *inputs = psArrayAlloc(argc); // Input filenames
     75    inputs->n = argc;
    7576    for (int i = 0; i < argc; i++) {
    76         inputs->data[i] = psAlloc(strlen(argv[i]));
    77         strncpy(inputs->data[i], argv[i], strlen(argv[i]));
     77        inputs->data[i] = psAlloc(strlen(argv[i]));
     78        strncpy(inputs->data[i], argv[i], strlen(argv[i]));
     79        psTrace("stac.size", 8, "Input file: %s\n", inputs->data[i]);
    7880    }
    7981
    8082    // Read input files
    81     psArray *images = stacReadImages(NULL, inputs);
     83    psVector *xSize = psVectorAlloc(inputs->n, PS_TYPE_S32);
     84    psVector *ySize = psVectorAlloc(inputs->n, PS_TYPE_S32);
     85    xSize->n = inputs->n;
     86    ySize->n = inputs->n;
     87    for (int i = 0; i < inputs->n; i++) {
     88        const char *filename = inputs->data[i]; // Name of FITS file
     89        psFits *file = psFitsOpen(filename, "r"); // FITS file handle
     90        psMetadata *header = psFitsReadHeader(NULL, file); // FITS header
     91        psFitsClose(file);
     92        bool mdok = true;               // Status of MD lookup
     93        xSize->data.S32[i] = psMetadataLookupS32(&mdok, header, "NAXIS1");
     94        if (!mdok) {
     95            psError(PS_ERR_IO, true, "Unable to find NAXIS1 in %s\n", filename);
     96            exit(EXIT_FAILURE);
     97        }
     98        ySize->data.S32[i] = psMetadataLookupS32(&mdok, header, "NAXIS2");
     99        if (!mdok) {
     100            psError(PS_ERR_IO, true, "Unable to find NAXIS2 in %s\n", filename);
     101            exit(EXIT_FAILURE);
     102        }
     103        psFree(header);
     104    }
    82105
    83106    // Read maps
    … …  
    85108
    86109    // Get size
    87     int outnx, outny;                   // Size of combined image
    88     float xMapDiff, yMapDiff;           // Difference to apply to maps
    89     stacSize(&outnx, &outny, &xMapDiff, &yMapDiff, images, maps);
     110    int outnx, outny;                   // Size of combined image
     111    float xMapDiff, yMapDiff;           // Difference to apply to maps
     112    stacSize(&outnx, &outny, &xMapDiff, &yMapDiff, xSize, ySize, maps);
     113    psFree(xSize);
     114    psFree(ySize);
    90115
    91116    printf("%d %d\n", outnx, outny);
    … …  
    95120
    96121    psFree(inputs);
    97     psFree(images);
    98122    psFree(maps);
    99123

  • trunk/stac/src/stac.c

    r6771 r6887  
    4545    // Generate masks
    4646    psArray *masks = psArrayAlloc(inputs->n);
     47    masks->n = inputs->n;
    4748    for (int i = 0; i < inputs->n; i++) {
    4849        psImage *image = inputs->data[i]; // Image for which to get mask
    … …  
    6768    // Get size, if not input
    6869    if (config->outnx == 0 || config->outny == 0) {
    69         stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps);
     70        psVector *xSize = psVectorAlloc(inputs->n, PS_TYPE_S32); // Sizes of images in x
     71        psVector *ySize = psVectorAlloc(inputs->n, PS_TYPE_S32); // Sizes of images in y
     72        xSize->n = inputs->n;
     73        ySize->n = inputs->n;
     74        for (int i = 0; i < inputs->n; i++) {
     75            psImage *image = inputs->data[i]; // The i-th image
     76            xSize->data.S32[i] = image->numCols;
     77            ySize->data.S32[i] = image->numRows;
     78        }
     79        stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, xSize, ySize, maps);
     80        psFree(xSize);
     81        psFree(ySize);
    7082    }
    7183
    … …  
    96108    for (int i = 0; i < inputs->n; i++) {
    97109        char errName[MAXCHAR];          // Filename of error image
    98         sprintf(errName,"%s.err",config->inputs->data[i]);
     110        sprintf(errName,"%s.err", (char*)config->inputs->data[i]);
    99111        psFits *errorFile = psFitsOpen(errName, "w");
    100         if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0)) {
     112        if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0, NULL)) {
    101113            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
    102114        }
    … …  
    118130        char shiftName[MAXCHAR];        // Filename of shift image
    119131        char errName[MAXCHAR];          // Filename of error image
    120         sprintf(shiftName,"%s.shift1",config->inputs->data[i]);
    121         sprintf(errName,"%s.shifterr1",config->inputs->data[i]);
     132        sprintf(shiftName,"%s.shift1", (char*)config->inputs->data[i]);
     133        sprintf(errName,"%s.shifterr1",(char*)config->inputs->data[i]);
    122134        psFits *shiftFile = psFitsOpen(shiftName, "w");
    123135        psFits *errFile = psFitsOpen(errName, "w");
    124136        psImage *trans = transformed->data[i]; // Transformed image
    125137        psImage *transErr = transformedErrors->data[i]; // Transformed error image
    126         if (!psFitsWriteImage(shiftFile, NULL, trans, 0)) {
     138        if (!psFitsWriteImage(shiftFile, NULL, trans, 0, NULL)) {
    127139            psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
    128140        }
    129141        psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
    130         if (!psFitsWriteImage(errFile, NULL, transErr, 0)) {
     142        if (!psFitsWriteImage(errFile, NULL, transErr, 0, NULL)) {
    131143            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
    132144        }
    … …  
    146158    psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits
    147159    psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits
     160    saturated->n = transformed->n;
     161    bad->n = transformed->n;
    148162    for (int i = 0; i < transformed->n; i++) {
    149163        saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
    … …  
    182196    for (int i = 0; i < rejected->n; i++) {
    183197        char rejName[MAXCHAR];  // Filename of rejection image
    184         sprintf(rejName, "%s.shiftrej", config->inputs->data[i]);
     198        sprintf(rejName, "%s.shiftrej", (char*)config->inputs->data[i]);
    185199
    186200        psFits *rejFile = psFitsOpen(rejName, "w");
    187         if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0)) {
     201        if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0, NULL)) {
    188202            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
    189203        }
    … …  
    196210    sprintf(preName, "%s.pre", config->output);
    197211    psFits *preFile = psFitsOpen(preName, "w");
    198     if (!psFitsWriteImage(preFile, NULL, combined, 0)) {
     212    if (!psFitsWriteImage(preFile, NULL, combined, 0, NULL)) {
    199213        psErrorStackPrint(stderr, "Unable to write image: %s\n", preName);
    200214    }
    … …  
    205219    // Get regions of interest in the source frame
    206220    psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest
     221    regions->n = inputs->n;
    207222    for (int i = 0; i < inputs->n; i++) {
    208223        regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
    … …  
    211226#ifdef TESTING
    212227        char regionName[MAXCHAR];       // Filename of region image
    213         sprintf(regionName,"%s.region",config->inputs->data[i]);
     228        sprintf(regionName,"%s.region",(char*)config->inputs->data[i]);
    214229        psFits *regionFile = psFitsOpen(regionName, "w");
    215         if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0)) {
     230        if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0, NULL)) {
    216231            psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName);
    217232        }

  • trunk/stac/src/stac.h

    r5743 r6887  
    1919// Read the input files and return an array of images
    2020psArray *stacReadImages(psArray **headers, // The image headers, to be returned
    21                         psArray *filenames // The file names of the images
     21                        psArray *filenames // The file names of the images
    2222    );
    2323
    … …  
    3838// Calculate the error images
    3939psArray *stacErrorImages(psArray *inputs, // Array of input images
    40                          float gain,    // Gain, in e/ADU
    41                          float rn       // Read noise, in e
     40                         float gain,    // Gain, in e/ADU
     41                         float rn       // Read noise, in e
    4242    );
    4343
    … …  
    4747
    4848// Transform input images, return success
    49 bool stacTransform(psArray **outputs,   // Transformed images for output
    50                    psArray **outErrors, // Transformed error images for output
    51                    const psArray *images, // Array of images to be transformed
    52                    const psArray *maps, // Array of polynomials that do the transformation
    53                    const psArray *errors, // Array of error images to be transformed
    54                    const psArray *masks, // Masks of input images
    55                    const psImage *region, // Region of interest for transformation
    56                    const psVector *scales, // Relative scales
    57                    const psVector *offsets, // Relative offsets
    58                    int outnx, int outny // Size of output images
     49bool stacTransform(psArray **outputs,   // Transformed images for output
     50                   psArray **outErrors, // Transformed error images for output
     51                   const psArray *images, // Array of images to be transformed
     52                   const psArray *maps, // Array of polynomials that do the transformation
     53                   const psArray *errors, // Array of error images to be transformed
     54                   const psArray *masks, // Masks of input images
     55                   const psImage *region, // Region of interest for transformation
     56                   const psVector *scales, // Relative scales
     57                   const psVector *offsets, // Relative offsets
     58                   int outnx, int outny // Size of output images
    5959    );
    6060
    … …  
    6868// Print out the problem
    6969void stacMemoryProblem(const psMemBlock* ptr, ///< the pointer to the problematic memory block.
    70                        const char *file, ///< the file in which the problem originated
    71                        psS32 lineno     ///< the line number in which the problem originated
     70                       const char *file, ///< the file in which the problem originated
     71                       psS32 lineno     ///< the line number in which the problem originated
    7272    );
    7373
    … …  
    8181
    8282// Get the mean for a bunch of values
    83 float stacCombineMean(psVector *values, // Values for which to take the mean
    84                       psVector *errors, // Errors in the values
    85                       psVector *masks   // Masks for the values, 0 = don't use, 1 = use
     83float stacCombineMean(psVector *values, // Values for which to take the mean
     84                      psVector *errors, // Errors in the values
     85                      psVector *masks   // Masks for the values, 0 = don't use, 1 = use
    8686    );
    8787
    8888// Get the median for a bunch of values
    8989float stacCombineMedian(psVector *values, // Values for which to take the median
    90                         psVector *errors, // Errors in the values
    91                         psVector *masks // Masks for the values, 0 = don't use, 1 = use
     90                        psVector *errors, // Errors in the values
     91                        psVector *masks // Masks for the values, 0 = don't use, 1 = use
    9292    );
    9393
    9494// Combine the transformed images
    95 bool stacCombine(psImage **combined,    // The combined image for output
    96                  psArray **rejected,    // Array of rejection masks
    97                  psArray *images,       // Array of transformed images
    98                  psArray *errors,       // Array of transformed error images
    99                  int nReject,           // Number of rejection iterations
    100                  psImage *region,       // Region to combine
    101                  psVector *saturated,   // Saturation limits for each image
    102                  psVector *bad,         // Bad pixel limits for each image
    103                  float reject           // Rejection (k-sigma)
     95bool stacCombine(psImage **combined,    // The combined image for output
     96                 psArray **rejected,    // Array of rejection masks
     97                 psArray *images,       // Array of transformed images
     98                 psArray *errors,       // Array of transformed error images
     99                 int nReject,           // Number of rejection iterations
     100                 psImage *region,       // Region to combine
     101                 psVector *saturated,   // Saturation limits for each image
     102                 psVector *bad,         // Bad pixel limits for each image
     103                 float reject           // Rejection (k-sigma)
    104104    );
    105105
    … …  
    110110// Invert an array of maps
    111111psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
    112                         int outnx, int outny // Size of output image
     112                        int outnx, int outny // Size of output image
    113113    );
    114114
    … …  
    118118
    119119// Return the relative gradient for a given pixel
    120 float stacGradient(psImage *image,      // Input for which to measure the gradient
    121                    int x, int y         // Coordinates at which to measure the gradient
     120float stacGradient(psImage *image,      // Input for which to measure the gradient
     121                   int x, int y         // Coordinates at which to measure the gradient
    122122    );
    123123
    124124// Transform the rejection masks back to the source frame
    125 psArray *stacRejection(psArray *inputs, // Input images
    126                        psArray *rejected, // Rejected images
    127                        psArray *regions, // Regions of interest
    128                        psArray *maps,   // Maps from input to transformed image
    129                        psArray *inverseMaps, // Maps from transformed to input image
    130                        float frac,      // Fraction of pixel rejected before looking more carefully
    131                        float grad,      // Gradient limit for rejection
    132                        psArray *names   // Names of original images (only for writing out when TESTING)
     125psArray *stacRejection(psArray *inputs, // Input images
     126                       psArray *rejected, // Rejected images
     127                       psArray *regions, // Regions of interest
     128                       psArray *maps,   // Maps from input to transformed image
     129                       psArray *inverseMaps, // Maps from transformed to input image
     130                       float frac,      // Fraction of pixel rejected before looking more carefully
     131                       float grad,      // Gradient limit for rejection
     132                       psArray *names   // Names of original images (only for writing out when TESTING)
    133133    );
    134134
    … …  
    139139// Returns the change in x' and y' for a change in x and y of 1
    140140bool stacPlaneTransformDeriv(psPlane *out, // Output derivative, assumed already allocated
    141                              psPlaneTransform *transform, // The transform for which to obtain the derivative
    142                              psPlane *in // The position at which to obtain the derivative
     141                             psPlaneTransform *transform, // The transform for which to obtain the derivative
     142                             psPlane *in // The position at which to obtain the derivative
    143143    );
    144144
    … …  
    146146// Basically, uses derivatives of the transformation to work out what pixels in the source might be of interest
    147147psImage *stacAreaOfInterest(psImage *mask, // Mask that is to be used to determine the area of interest
    148                             psPlaneTransform *map, // Map from the mask to the area of interest
    149                             int nxOut, int nyOut // Size of the area of interest
     148                            psPlaneTransform *map, // Map from the mask to the area of interest
     149                            int nxOut, int nyOut // Size of the area of interest
    150150    );
    151151
    … …  
    155155
    156156// Calculate the size of the output image
    157 bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
    158               float *xMapDiff, float *yMapDiff, // Difference applied to maps
    159               const psArray *images,    // Input images
    160               psArray *maps             // Transformation maps
     157bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
     158              float *xMapDiff, float *yMapDiff, // Difference applied to maps
     159              const psVector *xSizes,   // Sizes of images in x
     160              const psVector *ySizes,   // Sizes of images in y
     161              psArray *maps             // Transformation maps
    161162    );
    162163
    … …  
    167168// Calculate the background in an image
    168169float stacBackground(const psImage *image, // Image for which to get the background
    169                      int sample         // Sample in increments of this value
     170                     int sample         // Sample in increments of this value
    170171    );
    171172
    172173
    173174// Calculate the relative scales and offsets between the images
    174 bool stacScales(psVector **scalesPtr,   // Scales to return
    175                 psVector **offsetsPtr,  // Offsets to return
    176                 const psArray *images,  // Images on which to measure the scales and offsets
    177                 const char *starFile,   // File containing coordinates to photometer
    178                 const char *starMapFile, // Map for coodinates to the common output frame
    179                 float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
    180                 float aper              // Aperture to use for photometry (radius)
     175bool stacScales(psVector **scalesPtr,   // Scales to return
     176                psVector **offsetsPtr,  // Offsets to return
     177                const psArray *images,  // Images on which to measure the scales and offsets
     178                const char *starFile,   // File containing coordinates to photometer
     179                const char *starMapFile, // Map for coodinates to the common output frame
     180                float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
     181                float aper              // Aperture to use for photometry (radius)
    181182    );
    182183
    183184// Rescale images
    184 bool stacRescale(psArray *images,       // Images to rescale
    185                  psArray *errImages,    // Variance images to rescale
    186                  const psImage *mask,   // Mask indicating which pixels to scale
    187                 const psVector *scales,// Scales for images
    188                 const psVector *offsets // Offsets for images
     185bool stacRescale(psArray *images,       // Images to rescale
     186                 psArray *errImages,    // Variance images to rescale
     187                 const psImage *mask,   // Mask indicating which pixels to scale
     188                const psVector *scales,// Scales for images
     189                const psVector *offsets // Offsets for images
    189190    );
    190191
    … …  
    196197
    197198// Write map out
    198 bool stacWriteMap(const char *mapName,  // Filename to write to
    199                   psPlaneTransform *map // Map to write
     199bool stacWriteMap(const char *mapName,  // Filename to write to
     200                  psPlaneTransform *map // Map to write
    200201    );
    201202
    202203// Write multiple maps
    203204bool stacWriteMaps(const psArray *names, // Filenames of the input images (will add ".map")
    204                    const psArray *maps  // Maps to write
     205                   const psArray *maps  // Maps to write
    205206    );
    206207

  • trunk/stac/src/stacCombine.c

    r5745 r6887  
    103103    psVector *deltas = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the errors in the statistics step
    104104    psVector *mask = psVectorAlloc(nImages, PS_TYPE_U8); // Will hold the mask in the statistics step
     105    pixels->n = nImages;
     106    deltas->n = nImages;
     107    mask->n = nImages;
    105108
    106109    // Set up rejection masks
    … …  
    109112            // Allocate the rejection masks, if required
    110113            *rejected = psArrayAlloc(nImages);
     114            (*rejected)->n = nImages;
    111115        } else {
    112116            assert((*rejected)->n != nImages);
    … …  
    202206    char chiName[MAXCHAR];              // Filename of chi^2 image
    203207    sprintf(chiName,"chi2_%d.fits",iteration);
    204     psFits *chiFile = psFitsAlloc(chiName);
    205     if (!psFitsWriteImage(chiFile, NULL, chi2 , 0)) {
     208    psFits *chiFile = psFitsOpen(chiName, "w");
     209    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0, NULL)) {
    206210        psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
    207211    }
    208212    psTrace("stac.combine", 1, "Chi^2 image written to %s\n", chiName);
    209     psFree(chiFile);
     213    psFitsClose(chiFile);
    210214    psFree(chi2);
    211215#endif

  • trunk/stac/src/stacConfig.c

    r5743 r6887  
    2424    config->outnx = 0;
    2525    config->outny = 0;
    26     config->saturated = 65535.0;        // Saturation level
    27     config->bad = 0.0;                  // Bad level
     26    config->saturated = 65535.0;        // Saturation level
     27    config->bad = 0.0;                  // Bad level
    2828    config->reject = 2.5;
    2929    config->frac = 0.45;
    … …  
    3939    // Free the vectors, if necessary
    4040    if (config->inputs) {
    41         psFree(config->inputs);
     41        psFree(config->inputs);
    4242    }
    4343    // Free everything
    … …  
    4646
    4747
    48 stacConfig *stacParseConfig(int argc,   // Number of command-line arguments
    49                             char **argv // Command-line arguments
     48stacConfig *stacParseConfig(int argc,   // Number of command-line arguments
     49                            char **argv // Command-line arguments
    5050    )
    5151{
    5252    stacConfig *config = stacConfigAlloc(); // Configuration values
    53     const char *programName = argv[0];  // Program name
     53    const char *programName = argv[0];  // Program name
    5454
    5555    /* Variables for getopt */
    … …  
    6161    while ((opt = getopt(argc, argv, "hvSg:r:o:s:b:k:n:f:G:p:a:")) != -1) {
    6262        switch (opt) {
    63           case 'h':
    64             help(programName);
    65             exit(EXIT_SUCCESS);
    66           case 'v':
     63          case 'h':
     64            help(programName);
     65            exit(EXIT_SUCCESS);
     66          case 'v':
    6767            config->verbose++;
    6868            break;
    69           case 'g':
    70             if (sscanf(optarg, "%f", &config->gain) != 1) {
    71                 help(programName);
    72                 exit(EXIT_FAILURE);
    73             }
    74             break;
    75           case 'r':
    76             if (sscanf(optarg, "%f", &config->readnoise) != 1) {
    77                 help(programName);
    78                 exit(EXIT_FAILURE);
    79             }
    80             break;
    81           case 'o':
     69          case 'g':
     70            if (sscanf(optarg, "%f", &config->gain) != 1) {
     71                help(programName);
     72                exit(EXIT_FAILURE);
     73            }
     74            break;
     75          case 'r':
     76            if (sscanf(optarg, "%f", &config->readnoise) != 1) {
     77                help(programName);
     78                exit(EXIT_FAILURE);
     79            }
     80            break;
     81          case 'o':
    8282            if ((sscanf(argv[optind-1], "%d", &config->outnx) != 1) ||
    8383                (sscanf(argv[optind++], "%d", &config->outny) != 1)) {
    … …  
    8686                exit(EXIT_FAILURE);
    8787            }
    88             break;
    89           case 's':
    90             if (sscanf(optarg, "%f", &config->saturated) != 1) {
    91                 help(programName);
    92                 exit(EXIT_FAILURE);
    93             }
    94             break;
    95           case 'b':
    96             if (sscanf(optarg, "%f", &config->bad) != 1) {
    97                 help(programName);
    98                 exit(EXIT_FAILURE);
    99             }
    100             break;
    101           case 'p':
    102             if (argc < optind+1) {
    103                 help(programName);
    104                 exit(EXIT_FAILURE);
    105             }
    106             config->starFile = argv[optind-1];
    107             config->starMapFile = argv[optind++];
    108             // Note: incrementing optind, so I can read more than one parameter.
    109             break;
    110           case 'a':
    111             if (sscanf(optarg, "%f", &config->aper) != 1) {
    112                 help(programName);
    113                 exit(EXIT_FAILURE);
    114             }
    115             break;
    116           case 'k':
    117             if (sscanf(optarg, "%f", &config->reject) != 1) {
    118                 help(programName);
    119                 exit(EXIT_FAILURE);
    120             }
    121             break;
    122           case 'n':
    123             if (sscanf(optarg, "%d", &config->nReject) != 1) {
    124                 help(programName);
    125                 exit(EXIT_FAILURE);
    126             }
    127             break;
    128           case 'f':
    129             if (sscanf(optarg, "%f", &config->frac) != 1) {
    130                 help(programName);
    131                 exit(EXIT_FAILURE);
    132             }
    133             break;
    134           case 'G':
    135             if (sscanf(optarg, "%f", &config->grad) != 1) {
    136                 help(programName);
    137                 exit(EXIT_FAILURE);
    138             }
    139             break;
    140           case 'S':
    141             config->saveShifts = true;
    142             break;
    143           default:
    144             help(programName);
    145         }
     88            break;
     89          case 's':
     90            if (sscanf(optarg, "%f", &config->saturated) != 1) {
     91                help(programName);
     92                exit(EXIT_FAILURE);
     93            }
     94            break;
     95          case 'b':
     96            if (sscanf(optarg, "%f", &config->bad) != 1) {
     97                help(programName);
     98                exit(EXIT_FAILURE);
     99            }
     100            break;
     101          case 'p':
     102            if (argc < optind+1) {
     103                help(programName);
     104                exit(EXIT_FAILURE);
     105            }
     106            config->starFile = argv[optind-1];
     107            config->starMapFile = argv[optind++];
     108            // Note: incrementing optind, so I can read more than one parameter.
     109            break;
     110          case 'a':
     111            if (sscanf(optarg, "%f", &config->aper) != 1) {
     112                help(programName);
     113                exit(EXIT_FAILURE);
     114            }
     115            break;
     116          case 'k':
     117            if (sscanf(optarg, "%f", &config->reject) != 1) {
     118                help(programName);
     119                exit(EXIT_FAILURE);
     120            }
     121            break;
     122          case 'n':
     123            if (sscanf(optarg, "%d", &config->nReject) != 1) {
     124                help(programName);
     125                exit(EXIT_FAILURE);
     126            }
     127            break;
     128          case 'f':
     129            if (sscanf(optarg, "%f", &config->frac) != 1) {
     130                help(programName);
     131                exit(EXIT_FAILURE);
     132            }
     133            break;
     134          case 'G':
     135            if (sscanf(optarg, "%f", &config->grad) != 1) {
     136                help(programName);
     137                exit(EXIT_FAILURE);
     138            }
     139            break;
     140          case 'S':
     141            config->saveShifts = true;
     142            break;
     143          default:
     144            help(programName);
     145        }
    146146    }
    147147
    … …  
    154154    }
    155155
    156     config->output = argv[0];           // Output file
     156    config->output = argv[0];           // Output file
    157157    // Get the input files
    158     config->inputs = psArrayAlloc(argc-1);
     158    config->inputs = psArrayAlloc(argc - 1);
     159    config->inputs->n = argc - 1;
    159160    for (int i = 1; i < argc; i++) {
    160         config->inputs->data[i-1] = psAlloc(strlen(argv[i]));
    161         strncpy(config->inputs->data[i-1], argv[i], strlen(argv[i]));
     161        config->inputs->data[i-1] = psAlloc(strlen(argv[i]));
     162        strncpy(config->inputs->data[i-1], argv[i], strlen(argv[i]));
    162163    }
    163164
    … …  
    167168    psTrace("stac.config", 9, "%d inputs:\n",config->inputs->n);
    168169    for (int i = 0; i < config->inputs->n; i++) {
    169         psTrace("stac.config", 9, "\t%s\n", config->inputs->data[i]);
     170        psTrace("stac.config", 9, "\t%s\n", config->inputs->data[i]);
    170171    }
    171172    psTrace("stac.config", 9, "Output file is %s\n",config->output);

  • trunk/stac/src/stacErrorImages.c

    r3666 r6887  
    44
    55psArray *stacErrorImages(psArray *inputs, // Array of input images
    6                          float gain,    // Gain, in e/ADU
    7                          float rn       // Read noise, in e
     6                         float gain,    // Gain, in e/ADU
     7                         float rn       // Read noise, in e
    88    )
    99{
    10     float invGain = 1.0/gain;           // Inverse square root of gain
    11     rn /= gain;                         // Read noise in ADU
     10    float invGain = 1.0/gain;           // Inverse square root of gain
     11    rn /= gain;                         // Read noise in ADU
    1212    psArray *errors = psArrayAlloc(inputs->n);
     13    errors->n = inputs->n;
    1314
    1415    psTrace("stac.errors", 1, "Calculating error images....\n");
    … …  
    1617    // Iterate over the input images
    1718    for (int i = 0; i < inputs->n; i++) {
    18         psTrace("stac.errors",5,"Working on image #%d\n",i);
     19        psTrace("stac.errors",5,"Working on image #%d\n",i);
    1920
    20         psImage *image = inputs->data[i]; // Pull out the image of interest
    21         int numRows = image->numRows;   // Number of rows
    22         int numCols = image->numCols;   // Number of columns   
    23         psImage *error = psImageAlloc(numCols, numRows, PS_TYPE_F32); // The error image
     21        psImage *image = inputs->data[i]; // Pull out the image of interest
     22        int numRows = image->numRows;   // Number of rows
     23        int numCols = image->numCols;   // Number of columns
     24        psImage *error = psImageAlloc(numCols, numRows, PS_TYPE_F32); // The error image
    2425
    25         // Iterate over the pixels
    26         for (int r = 0; r < numRows; r++) {
    27             for (int c = 0; c < numCols; c++) {
    28                 // We actually calculate the variance
    29                 error->data.F32[r][c] = image->data.F32[r][c]*invGain + rn*rn;
    30             }
    31         }
     26        // Iterate over the pixels
     27        for (int r = 0; r < numRows; r++) {
     28            for (int c = 0; c < numCols; c++) {
     29                // We actually calculate the variance
     30                error->data.F32[r][c] = image->data.F32[r][c]*invGain + rn*rn;
     31            }
     32        }
    3233
    33         // Put image onto the array
    34         errors->data[i] = error;
     34        // Put image onto the array
     35        errors->data[i] = error;
    3536    }
    3637

  • trunk/stac/src/stacInvertMaps.c

    r6450 r6887  
    1414    int nMaps = maps->n;                // Number of maps
    1515    psArray *inverted = psArrayAlloc(nMaps); // Array of inverted maps for output
     16    inverted->n = nMaps;
    1617
    1718    psTrace("stac.invertMaps", 1, "Inverting maps....\n");
    … …  
    4546        psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
    4647        psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     48        xIn->n = NUM_GRID * NUM_GRID;
     49        yIn->n = NUM_GRID * NUM_GRID;
     50        xOut->n = NUM_GRID * NUM_GRID;
     51        yOut->n = NUM_GRID * NUM_GRID;
    4752
    4853        // Create grid of points
    … …  
    6671        psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
    6772        psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
     73        xVector->n = nCoeff;
     74        yVector->n = nCoeff;
    6875        for (int i = 0; i < nCoeff; i++) {
    6976            for (int j = 0; j < nCoeff; j++) {
    … …  
    104111        // Solution via LU Decomposition
    105112        psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     113        permutation->n = nCoeff;
    106114        psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
    107115        psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x

  • trunk/stac/src/stacRead.c

    r6771 r6887  
    1313    int nFiles = filenames->n;          // The number of input files
    1414    psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
     15    images->n = nFiles;
    1516    assert(!headers || ! *headers || (*headers)->n == nFiles);
    1617    if (headers && ! *headers) {
    1718        *headers = psArrayAlloc(nFiles);
     19        (*headers)->n = nFiles;
    1820    }
    1921
    … …  
    7173    psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
    7274    float x, y;                         // Coordinates to read
    73     int num = 0;                        // Number of coordinates read
    7475    while (fscanf(file, "%f %f\n", &x, &y) == 2) {
    7576        psPlane *coord = psPlaneAlloc();// A coordinate
    7677        coord->x = x;
    7778        coord->y = y;
    78         coords->data[num] = coord;
    79         num++;
    80         if (num % BUFFER) {
    81             coords = psArrayRealloc(coords, num + BUFFER);
    82         }
    83     }
    84     coords->n = num;
     79        coords->data[coords->n] = coord;
     80        coords->n++;
     81        if (coords->n % BUFFER == 0) {
     82            coords = psArrayRealloc(coords, coords->n + BUFFER);
     83        }
     84    }
    8585
    8686    psTrace("stac.read.coords", 5, "%d coordinates read.\n", num);
    … …  
    184184    int nFiles = filenames->n;          // The number of input files
    185185    psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
     186    maps->n = nFiles;
    186187    char mapfile[MAXCHAR];              // Filename of map
    187188

  • trunk/stac/src/stacRejection.c

    r5743 r6887  
    88#define MIN(x,y) ((x) < (y) ? (x) : (y))
    99
    10 float stacGradient(psImage *image,      // Input for which to measure the gradient
    11                    int x, int y         // Coordinates at which to measure the gradient
     10float stacGradient(psImage *image,      // Input for which to measure the gradient
     11                   int x, int y         // Coordinates at which to measure the gradient
    1212    )
    1313{
    … …  
    2222    int yMax = MIN(y + 1, image->numRows - 1);
    2323    for (int j = yMin; j <= yMax; j++) {
    24         for (int i = xMin; i <= xMax; i++) {
    25             if ((i != x) && (j != y)) {
    26                 pixels->data.F32[num] = image->data.F32[j][i];
    27                 mask->data.U8[num] = 0;
    28                 num++;
    29             }
    30         }
     24        for (int i = xMin; i <= xMax; i++) {
     25            if ((i != x) && (j != y)) {
     26                pixels->data.F32[num] = image->data.F32[j][i];
     27                mask->data.U8[num] = 0;
     28                num++;
     29            }
     30        }
    3131    }
    3232    pixels->n = num;
    … …  
    4343}
    4444
    45 psArray *stacRejection(psArray *inputs, // Input images
    46                        psArray *rejected, // Rejected images
    47                        psArray *regions, // Regions of interest
    48                        psArray *maps,   // Maps from input to transformed image
    49                        psArray *inverseMaps, // Maps from transformed to input image
    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)
     45psArray *stacRejection(psArray *inputs, // Input images
     46                       psArray *rejected, // Rejected images
     47                       psArray *regions, // Regions of interest
     48                       psArray *maps,   // Maps from input to transformed image
     49                       psArray *inverseMaps, // Maps from transformed to input image
     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)
    5353    )
    5454{
    55     int nImages = inputs->n;            // Number of input images
     55    int nImages = inputs->n;            // Number of input images
    5656
    5757    // Check inputs
    … …  
    6363
    6464    for (int i = 0; i < nImages; i++) {
    65         psImage *input = inputs->data[i];
    66         psImage *region = regions->data[i];
    67         assert(input->numRows == region->numRows && input->numCols == region->numCols);
     65        psImage *input = inputs->data[i];
     66        psImage *region = regions->data[i];
     67        assert(input->numRows == region->numRows && input->numCols == region->numCols);
    6868    }
    6969
    … …  
    7777    psVector *grads = psVectorAlloc(nImages, PS_TYPE_F32); // Gradient for each image
    7878    psVector *gradsMask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask for gradient vector
     79    grads->n = nImages;
     80    gradsMask->n = nImages;
    7981
    8082    // Transform rejection masks back to source
    8183    psArray *inputRej = psArrayAlloc(nImages);
     84    inputRej->n = nImages;
    8285    for (int i = 0; i < nImages; i++) {
    83         // Size of input image
    84         int nxInput = ((psImage*)(inputs->data[i]))->numCols;
    85         int nyInput = ((psImage*)(inputs->data[i]))->numRows;
    86         psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_U8); // The pixel mask for the input
    87 #ifdef TESTING
    88         psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
    89         psImage *gradient = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);        // The gradient image
    90 #endif
    91         psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
    92         psPlaneTransform *map = maps->data[i]; // The map from input to output
    93         psImage *region = regions->data[i]; // The region of interest for this image
    94 
    95         psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
    96         int nBad = 0;                   // Number of bad pixels
     86        // Size of input image
     87        int nxInput = ((psImage*)(inputs->data[i]))->numCols;
     88        int nyInput = ((psImage*)(inputs->data[i]))->numRows;
     89        psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_U8); // The pixel mask for the input
     90#ifdef TESTING
     91        psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
     92        psImage *gradient = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);        // The gradient image
     93#endif
     94        psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
     95        psPlaneTransform *map = maps->data[i]; // The map from input to output
     96        psImage *region = regions->data[i]; // The region of interest for this image
     97
     98        psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
     99        int nBad = 0;                   // Number of bad pixels
    97100
    98101#ifdef CRFLUX
    99         // Set up CR output
    100         FILE *crs = NULL;               // File for outputting details of rejected pixels
    101         char crfile[MAXCHAR];   // Filename
    102         sprintf(crfile,"%s.cr",names->data[i]);
    103         if ((crs = fopen(crfile, "w")) == NULL) {
    104             fprintf(stderr, "Unable to open file for detailed output, %s\n");
    105             return NULL;
    106         }
    107 #endif
    108 
    109         // Transform the mask
    110         // Optimisation option is to only transform the pixels that have been rejected in the output,
    111         // calculate derivatives of the map, and use that as a buffer around the transformed position
    112         // in the input image.
    113         psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    114         for (int y = 0; y < nyInput; y++) {
    115             for (int x = 0; x < nxInput; x++) {
    116 
    117                 // Only transform pixels of interest
    118                 if (region->data.U8[y][x]) {
    119 
    120                     inCoords->x = (double)x + 0.5;
    121                     inCoords->y = (double)y + 0.5;
    122                     (void)psPlaneTransformApply(outCoords, map, inCoords);
    123                     float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
    124                                                                    NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
    125 #ifdef TESTING
    126                     rejmap->data.F32[y][x] = maskVal;
    127 #endif
    128 
    129                     if (maskVal > frac) {
    130                         // Calculate mean gradient on other images
    131                         float meanGrads = 0.0;
    132                         int numGrads = 0;
    133                         for (int j = 0; j < nImages; j++) {
    134                             if (i != j) {
    135                                 // Get coordinates for that image
    136                                 (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
    137                                 int xPix = (int)(inCoords->x + 0.5);
    138                                 int yPix = (int)(inCoords->y + 0.5);
    139                                 if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
    140                                     (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
    141                                     // Calculate the gradient
    142                                     grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
    143                                     gradsMask->data.U8[j] = 0;
    144                                     numGrads++;
    145                                 } else {
    146                                     gradsMask->data.U8[j] = 1; // Mask this one
    147                                 }
    148                             } else {
    149                                 gradsMask->data.U8[j] = 1; // Mask this one
    150                             }
    151                         }
    152                         if (numGrads > 0) {
    153                             (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
    154                             meanGrads = median->sampleMedian;
    155                         } else {
    156                             meanGrads = 0.0;
    157                         }
    158 
    159 #ifdef TESTING
    160                         //gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
    161                         gradient->data.F32[y][x] = meanGrads;
    162 #endif
    163 
    164                         //if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
    165                         if (meanGrads > grad) {
    166                             mask->data.U8[y][x] = 1;
    167                             nBad++;
     102        // Set up CR output
     103        FILE *crs = NULL;               // File for outputting details of rejected pixels
     104        char crfile[MAXCHAR];   // Filename
     105        sprintf(crfile,"%s.cr",names->data[i]);
     106        if ((crs = fopen(crfile, "w")) == NULL) {
     107            fprintf(stderr, "Unable to open file for detailed output, %s\n");
     108            return NULL;
     109        }
     110#endif
     111
     112        // Transform the mask
     113        // Optimisation option is to only transform the pixels that have been rejected in the output,
     114        // calculate derivatives of the map, and use that as a buffer around the transformed position
     115        // in the input image.
     116        psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     117        for (int y = 0; y < nyInput; y++) {
     118            for (int x = 0; x < nxInput; x++) {
     119
     120                // Only transform pixels of interest
     121                if (region->data.U8[y][x]) {
     122
     123                    inCoords->x = (double)x + 0.5;
     124                    inCoords->y = (double)y + 0.5;
     125                    (void)psPlaneTransformApply(outCoords, map, inCoords);
     126                    float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
     127                                                                   NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     128#ifdef TESTING
     129                    rejmap->data.F32[y][x] = maskVal;
     130#endif
     131
     132                    if (maskVal > frac) {
     133                        // Calculate mean gradient on other images
     134                        float meanGrads = 0.0;
     135                        int numGrads = 0;
     136                        for (int j = 0; j < nImages; j++) {
     137                            if (i != j) {
     138                                // Get coordinates for that image
     139                                (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
     140                                int xPix = (int)(inCoords->x + 0.5);
     141                                int yPix = (int)(inCoords->y + 0.5);
     142                                if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
     143                                    (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
     144                                    // Calculate the gradient
     145                                    grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
     146                                    gradsMask->data.U8[j] = 0;
     147                                    numGrads++;
     148                                } else {
     149                                    gradsMask->data.U8[j] = 1; // Mask this one
     150                                }
     151                            } else {
     152                                gradsMask->data.U8[j] = 1; // Mask this one
     153                            }
     154                        }
     155                        if (numGrads > 0) {
     156                            (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
     157                            meanGrads = median->sampleMedian;
     158                        } else {
     159                            meanGrads = 0.0;
     160                        }
     161
     162#ifdef TESTING
     163                        //gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
     164                        gradient->data.F32[y][x] = meanGrads;
     165#endif
     166
     167                        //if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
     168                        if (meanGrads > grad) {
     169                            mask->data.U8[y][x] = 1;
     170                            nBad++;
    168171
    169172#ifdef CRFLUX
    170                             fprintf(crs, "%d %d --> %f %f %f\n", x, y,
    171                                     ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal,
    172                                     stacGradient(inputs->data[i], x, y));
    173 #endif
    174 
    175                         } else {
    176                             mask->data.U8[y][x] = 0;
    177                         } // Gradient threshold
    178                     } else {
    179                         mask->data.U8[y][x] = 0;
    180                     } // Rejection threshold
    181 
    182                 } else {
    183                     mask->data.U8[y][x] = 0;
    184                 } // Only touching pixels of interest
    185 
    186             }
    187         } // Iterating over pixels
    188         psFree(median);
     173                            fprintf(crs, "%d %d --> %f %f %f\n", x, y,
     174                                    ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal,
     175                                    stacGradient(inputs->data[i], x, y));
     176#endif
     177
     178                        } else {
     179                            mask->data.U8[y][x] = 0;
     180                        } // Gradient threshold
     181                    } else {
     182                        mask->data.U8[y][x] = 0;
     183                    } // Rejection threshold
     184
     185                } else {
     186                    mask->data.U8[y][x] = 0;
     187                } // Only touching pixels of interest
     188
     189            }
     190        } // Iterating over pixels
     191        psFree(median);
    189192
    190193#ifdef CRFLUX
    191         // Close file used for detailed output
    192         fclose(crs);
    193 #endif
    194 
    195         psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
    196         // Clip the image, and convert to suitable mask format
    197 
    198 #ifdef TESTING
    199         // Write error image out to check
    200         char maskName[MAXCHAR];         // Filename of mask image
    201         char rejmapName[MAXCHAR];       // Filename of rejection image
    202         char gradName[MAXCHAR];         // Filename of gradient image
    203         sprintf(maskName, "%s.mask", names->data[i]);
    204         sprintf(rejmapName, "%s.rejmap", names->data[i]);
    205         sprintf(gradName, "%s.grad", names->data[i]);
    206 
    207         psFits *maskFile = psFitsAlloc(maskName);
    208         psFits *rejmapFile = psFitsAlloc(rejmapName);
    209         psFits *gradFile = psFitsAlloc(gradName);
    210         if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
    211             psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
    212         }
    213         psTrace("stac", 1, "Mask image written to %s\n", maskName);
    214         if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0)) {
    215             psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
    216         }
    217         psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
    218         if (!psFitsWriteImage(gradFile, NULL, gradient, 0)) {
    219             psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
    220         }
    221         psTrace("stac", 1, "Gradient image written to %s\n", gradName);
    222         psFree(maskFile);
    223         psFree(rejmapFile);
    224         psFree(gradFile);
    225         psFree(rejmap);
    226         psFree(gradient);
    227 #endif
    228 
    229         // Stuff into the array
    230         inputRej->data[i] = mask;
     194        // Close file used for detailed output
     195        fclose(crs);
     196#endif
     197
     198        psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
     199        // Clip the image, and convert to suitable mask format
     200
     201#ifdef TESTING
     202        // Write error image out to check
     203        char maskName[MAXCHAR];         // Filename of mask image
     204        char rejmapName[MAXCHAR];       // Filename of rejection image
     205        char gradName[MAXCHAR];         // Filename of gradient image
     206        sprintf(maskName, "%s.mask", (char*)names->data[i]);
     207        sprintf(rejmapName, "%s.rejmap", (char*)names->data[i]);
     208        sprintf(gradName, "%s.grad", (char*)names->data[i]);
     209
     210        psFits *maskFile = psFitsOpen(maskName, "w");
     211        psFits *rejmapFile = psFitsOpen(rejmapName, "w");
     212        psFits *gradFile = psFitsOpen(gradName, "w");
     213        if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
     214            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
     215        }
     216        psTrace("stac", 1, "Mask image written to %s\n", maskName);
     217        if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0, NULL)) {
     218            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
     219        }
     220        psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
     221        if (!psFitsWriteImage(gradFile, NULL, gradient, 0, NULL)) {
     222            psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
     223        }
     224        psTrace("stac", 1, "Gradient image written to %s\n", gradName);
     225        psFitsClose(maskFile);
     226        psFitsClose(rejmapFile);
     227        psFitsClose(gradFile);
     228        psFree(rejmap);
     229        psFree(gradient);
     230#endif
     231
     232        // Stuff into the array
     233        inputRej->data[i] = mask;
    231234
    232235    }

  • trunk/stac/src/stacScales.c

    r5745 r6887  
    1919    psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample
    2020    psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
     21    values->n = numSamples + 1;
     22    mask->n = numSamples + 1;
    2123
    2224    int offset = 0;                     // Offset from start of the row
    … …  
    8486    } else {
    8587        scales = psVectorAlloc(images->n, PS_TYPE_F32);
     88        scales->n = images->n;
    8689        *scalesPtr = scales;
    8790    }
    … …  
    9396    } else {
    9497        offsets = psVectorAlloc(images->n, PS_TYPE_F32);
     98        offsets->n = images->n;
    9599        *offsetsPtr = offsets;
    96100    }
    … …  
    122126        // Transform the stellar positions to match the transformed reference frame
    123127        psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
     128        starCoordsTransformed->n = starCoords->n;
    124129        // Fix up difference between map and output frame
    125130        starMap->x->coeff[0][0] -= xMapDiff;
    … …  
    132137        psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
    133138        psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
     139        stars->n = images->n;
     140        masks->n = images->n;
    134141
    135142        // Set scales relative to the first image
    … …  
    144151            psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
    145152            psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
     153            photometry->n = starCoords->n;
     154            mask->n = starCoords->n;
    146155            for (int j = 0; j < starCoordsTransformed->n; j++) {
    147156                psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star

  • trunk/stac/src/stacSize.c

    r3683 r6887  
    44#include "stac.h"
    55
    6 bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
    7               float *xMapDiff, float *yMapDiff, // Difference applied to maps
    8               const psArray *images,    // Input images
    9               psArray *maps             // Transformation maps
     6bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
     7              float *xMapDiff, float *yMapDiff, // Difference applied to maps
     8              const psVector *xSizes,   // Sizes of images in x
     9              const psVector *ySizes,   // Sizes of images in y
     10              psArray *maps             // Transformation maps
    1011    )
    1112{
    12     int nImages = images->n;            // Number of input images
     13    int nImages = xSizes->n;            // Number of input images
     14    assert(ySizes->n == nImages);
     15    assert(xSizes->type.type == PS_TYPE_S32);
     16    assert(ySizes->type.type == PS_TYPE_S32);
    1317    assert(maps->n == nImages);
    14     assert(outnx && outny);             // Must be able to write to these
     18    assert(outnx && outny);             // Must be able to write to these
    1519
    1620    psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates
    … …  
    2428
    2529    for (int i = 0; i < nImages; i++) {
    26         psImage *image = images->data[i]; // The image
    27         psPlaneTransform *map = maps->data[i]; // The map
     30        int xSize = xSizes->data.S32[i];// Size of image in x
     31        int ySize = ySizes->data.S32[i];// Size of image in y
     32        psPlaneTransform *map = maps->data[i]; // The map
    2833
    29         psTrace("stac.size", 4, "Image %d:\n", i);
     34        psTrace("stac.size", 4, "Image %d:\n", i);
    3035
    31         // Lower left corner
    32         inCoord->x = 0.0;
    33         inCoord->y = 0.0;
    34         (void)psPlaneTransformApply(outCoord, map, inCoord);
    35         if (outCoord->x < xMin) {
    36             xMin = outCoord->x;
    37         } else if (outCoord->x > xMax) {
    38             xMax = outCoord->x;
    39         }
    40         if (outCoord->y < yMin) {
    41             yMin = outCoord->y;
    42         } else if (outCoord->y > yMax) {
    43             yMax = outCoord->y;
    44         }
    45         psTrace("stac.size", 4, "Lower left: %f %f\n", outCoord->x, outCoord->y);
     36        // Lower left corner
     37        inCoord->x = 0.0;
     38        inCoord->y = 0.0;
     39        (void)psPlaneTransformApply(outCoord, map, inCoord);
     40        if (outCoord->x < xMin) {
     41            xMin = outCoord->x;
     42        } else if (outCoord->x > xMax) {
     43            xMax = outCoord->x;
     44        }
     45        if (outCoord->y < yMin) {
     46            yMin = outCoord->y;
     47        } else if (outCoord->y > yMax) {
     48            yMax = outCoord->y;
     49        }
     50        psTrace("stac.size", 4, "Lower left: %f %f\n", outCoord->x, outCoord->y);
    4651
    47         // Lower right corner
    48         inCoord->x = image->numCols;
    49         (void)psPlaneTransformApply(outCoord, map, inCoord);
    50         if (outCoord->x < xMin) {
    51             xMin = outCoord->x;
    52         } else if (outCoord->x > xMax) {
    53             xMax = outCoord->x;
    54         }
    55         if (outCoord->y < yMin) {
    56             yMin = outCoord->y;
    57         } else if (outCoord->y > yMax) {
    58             yMax = outCoord->y;
    59         }
    60         psTrace("stac.size", 4, "Lower right: %f %f\n", outCoord->x, outCoord->y);
     52        // Lower right corner
     53        inCoord->x = xSize;
     54        (void)psPlaneTransformApply(outCoord, map, inCoord);
     55        if (outCoord->x < xMin) {
     56            xMin = outCoord->x;
     57        } else if (outCoord->x > xMax) {
     58            xMax = outCoord->x;
     59        }
     60        if (outCoord->y < yMin) {
     61            yMin = outCoord->y;
     62        } else if (outCoord->y > yMax) {
     63            yMax = outCoord->y;
     64        }
     65        psTrace("stac.size", 4, "Lower right: %f %f\n", outCoord->x, outCoord->y);
    6166
    62         // Upper right corner
    63         inCoord->x = 0.0;
    64         inCoord->y = image->numRows;
    65         (void)psPlaneTransformApply(outCoord, map, inCoord);
    66         if (outCoord->x < xMin) {
    67             xMin = outCoord->x;
    68         } else if (outCoord->x > xMax) {
    69             xMax = outCoord->x;
    70         }
    71         if (outCoord->y < yMin) {
    72             yMin = outCoord->y;
    73         } else if (outCoord->y > yMax) {
    74             yMax = outCoord->y;
    75         }
    76         psTrace("stac.size", 4, "Upper right: %f %f\n", outCoord->x, outCoord->y);
     67        // Upper right corner
     68        inCoord->x = 0.0;
     69        inCoord->y = ySize;
     70        (void)psPlaneTransformApply(outCoord, map, inCoord);
     71        if (outCoord->x < xMin) {
     72            xMin = outCoord->x;
     73        } else if (outCoord->x > xMax) {
     74            xMax = outCoord->x;
     75        }
     76        if (outCoord->y < yMin) {
     77            yMin = outCoord->y;
     78        } else if (outCoord->y > yMax) {
     79            yMax = outCoord->y;
     80        }
     81        psTrace("stac.size", 4, "Upper right: %f %f\n", outCoord->x, outCoord->y);
    7782
    78         // Upper left corner
    79         inCoord->x = image->numCols;
    80         (void)psPlaneTransformApply(outCoord, map, inCoord);
    81         if (outCoord->x < xMin) {
    82             xMin = outCoord->x;
    83         } else if (outCoord->x > xMax) {
    84             xMax = outCoord->x;
    85         }
    86         if (outCoord->y < yMin) {
    87             yMin = outCoord->y;
    88         } else if (outCoord->y > yMax) {
    89             yMax = outCoord->y;
    90         }
    91         psTrace("stac.size", 4, "Upper left: %f %f\n", outCoord->x, outCoord->y);
     83        // Upper left corner
     84        inCoord->x = xSize;
     85        (void)psPlaneTransformApply(outCoord, map, inCoord);
     86        if (outCoord->x < xMin) {
     87            xMin = outCoord->x;
     88        } else if (outCoord->x > xMax) {
     89            xMax = outCoord->x;
     90        }
     91        if (outCoord->y < yMin) {
     92            yMin = outCoord->y;
     93        } else if (outCoord->y > yMax) {
     94            yMax = outCoord->y;
     95        }
     96        psTrace("stac.size", 4, "Upper left: %f %f\n", outCoord->x, outCoord->y);
    9297    }
    9398
    9499    // Tweak the maps to account for the offset
    95100    if (xMapDiff != NULL) {
    96         *xMapDiff = floorf(xMin);
     101        *xMapDiff = floorf(xMin);
    97102    }
    98103    if (yMapDiff != NULL) {
    99         *yMapDiff = floorf(yMin);
     104        *yMapDiff = floorf(yMin);
    100105    }
    101106    for (int i = 0; i < nImages; i++) {
    102         psPlaneTransform *map = maps->data[i]; // The map of interest
    103         map->x->coeff[0][0] -= floorf(xMin);
    104         map->y->coeff[0][0] -= floorf(yMin);
     107        psPlaneTransform *map = maps->data[i]; // The map of interest
     108        map->x->coeff[0][0] -= floorf(xMin);
     109        map->y->coeff[0][0] -= floorf(yMin);
    105110    }
    106111

  • trunk/stac/src/stacTransform.c

    r5745 r6887  
    127127    if (*outputs == NULL) {
    128128        *outputs = psArrayAlloc(nImages);
     129        (*outputs)->n = nImages;
    129130        psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
    130131        for (int i = 0; i < nImages; i++) {
    … …  
    138139    if (errors && (*outErrors == NULL)) {
    139140        *outErrors = psArrayAlloc(errors->n);
     141        (*outErrors)->n = errors->n;
    140142        psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
    141143        for (int i = 0; i < nImages; i++) {

  • trunk/stac/src/sum.c

    r6771 r6887  
    1010    const char *outputName = argv[1];   // Output file name
    1111    psArray *inputNames = psArrayAlloc(argc-2);
     12    inputNames->n = argc - 2;
    1213    for (int i = 2; i < argc; i++) {
    1314        inputNames->data[i-2] = psStringCopy(argv[i]);
Note: See TracChangeset for help on using the changeset viewer.