IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5743


Ignore:
Timestamp:
Dec 7, 2005, 2:57:14 PM (20 years ago)
Author:
Paul Price
Message:

Importing from PAP cvs tree again. This will now be the working tree.

Location:
trunk/stac/src
Files:
3 added
1 deleted
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/Makefile

    r3683 r5743  
    11SHELL = /bin/sh
    22CC = gcc
    3 CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -D_GNU_SOURCE # -DTESTING -DCRFLUX
    4 PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2 -lmysqlclient
    5 LDFLAGS += $(PSLIB)
    63
     4### psLib
     5PSLIB_INCLUDE := $(shell pslib-config --cflags)
     6PSLIB_LINK := $(shell pslib-config --libs)
     7
     8### Additional flags for diagnostics, testing, etc
     9#PAP_CFLAGS = -DTESTING
     10#PAP_CFLAGS = -DCRFLUX
     11
     12### Build flags
     13CFLAGS += -g -O2 -std=c99 -Werror -D_GNU_SOURCE -DPS_NO_TRACE $(PSLIB_INCLUDE) $(PAP_CFLAGS)
     14LDFLAGS += $(PSLIB_LINK)
     15
     16### Ingredients for each program
    717STAC = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o stacHelp.o \
    8         stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o stacScales.o time.o
     18        stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o stacScales.o \
     19        stacTime.o
    920
    1021SHIFT = shift.o stacRead.o stacTransform.o stacCheckMemory.o stacInvertMaps.o stacSize.o
    1122
    1223COMBINE = combine.o combineConfig.o stacRead.o stacErrorImages.o stacCombine.o stacScales.o \
    13         stacCheckMemory.o time.o
     24        stacCheckMemory.o stacTime.o
    1425
    1526SHIFTSIZE = shiftSize.o stacWrite.o stacConfig.o stacRead.o stacCheckMemory.o stacSize.o
    … …  
    1728FIXIMAGE = fixImage.o
    1829
    19 TARGETS = stac shift combine shiftSize fixImage
     30CALCGRADIENT = calcGradient.o stacRejection.o
    2031
     32SUM = sum.o
     33
     34### List of targets
     35TARGETS = stac shift combine shiftSize calcGradient sum
     36
     37
     38### Build recipes
    2139.PHONY: tags clean empty test profile optimise all
    2240
    2341.c.o:
    2442                $(CC) -c $(CFLAGS) $(OPTFLAGS) $<
     43
     44all:            $(TARGETS)
    2545
    2646stac:           $(STAC)
    … …  
    3656                $(CC) $(CFLAGS) -o $@ $(SHIFTSIZE) $(LDFLAGS) $(OPTFLAGS)
    3757
    38 fixImage:       $(FIXIMAGE)
    39                 $(CC) $(CFLAGS) -o $@ $(FIXIMAGE) $(LDFLAGS) $(OPTFLAGS)
     58calcGradient:   $(CALCGRADIENT)
     59                $(CC) $(CFLAGS) -o $@ $(CALCGRADIENT) $(LDFLAGS) $(OPTFLAGS)
    4060
    41 all:            $(TARGETS)
     61sum:            $(SUM)
     62                $(CC) $(CFLAGS) -o $@ $(SUM) $(LDFLAGS) $(OPTFLAGS)
     63
    4264
    4365clean:

  • trunk/stac/src/combine.c

    r3681 r5743  
    2121    combineConfig *config = combineParseConfig(argc, argv); // Configuration
    2222
    23     psArray *images = stacReadImages(config->inNames); // The images
     23    psArray *headers = NULL;            // Array of headers
     24    psArray *images = stacReadImages(&headers, config->inNames); // The images
    2425    // Make fake errors
    2526    psArray *errors = psArrayAlloc(images->n);
    2627    for (int i = 0; i < images->n; i++) {
    27         psImage *image = images->data[i];
    28         psImage *err = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
    29         for (int y = 0; y < image->numRows; y++) {
    30             for (int x = 0; x < image->numCols; x++) {
    31                 err->data.F32[y][x] = 1.0;
    32             }
    33         }
    34         errors->data[i] = err;
     28        psImage *image = images->data[i];
     29        psImage *err = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     30        for (int y = 0; y < image->numRows; y++) {
     31            for (int x = 0; x < image->numCols; x++) {
     32                err->data.F32[y][x] = 1.0;
     33            }
     34        }
     35        errors->data[i] = err;
    3536    }
    3637
    3738    // Calculate scales between images
    38     psVector *scales = NULL;            // Relative scales between images
    39     psVector *offsets = NULL;           // Offsets between images
     39    psVector *scales = NULL;            // Relative scales between images
     40    psVector *offsets = NULL;           // Offsets between images
    4041    (void)stacScales(&scales, &offsets, images, config->starFile, config->starMap, 0.0, 0.0, config->aper);
    4142
    … …  
    4445    psVector *bad = psVectorAlloc(images->n, PS_TYPE_F32); // Bad limits
    4546    for (int i = 0; i < images->n; i++) {
    46         saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
    47         bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
     47        saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
     48        bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
    4849    }
    4950
    … …  
    5253
    5354    // Combine the images
    54     psImage *combined = NULL;           // Combined image
    55     psArray *rejected = NULL;           // Array of rejection masks
     55    psImage *combined = NULL;           // Combined image
     56    psArray *rejected = NULL;           // Array of rejection masks
    5657    stacCombine(&combined, &rejected, images, errors, config->nReject, NULL, saturated, bad, config->reject);
    5758
    58     psFits *outFile = psFitsAlloc(config->outName);
    59     if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) {
    60         psErrorStackPrint(stderr, "Unable to write image: %s\n", config->outName);
     59    psFits *outFile = psFitsOpen(config->outName, "r");
     60    if (!psFitsWriteImage(outFile, headers->data[0], combined, 0)) {
     61        psErrorStackPrint(stderr, "Unable to write image: %s\n", config->outName);
    6162    }
    6263    psTrace("combine", 1, "Combined image written to %s\n", config->outName);
    63     psFree(outFile);
     64    psFitsClose(outFile);
    6465
    6566    // Clean up
     67    psFree(headers);
    6668    psFree(combined);
    6769    psFree(rejected);

  • trunk/stac/src/shift.c

    r3673 r5743  
    1111{
    1212    fprintf (stderr, "shift: shift an image, given the transformation\n"
    13              "Usage: %s [-h] [-v] [-s NX NY] IN OUT\n"
    14              "where\n"
    15              "\t-h           Help (this info)\n"
    16              "\t-v           Increase verbosity level\n"
    17              "\t-s NX NY     Size of output image\n"
    18              "\tIN           Input image, which has associated .map file\n"
    19              "\tOUT          Output image\n",
    20              programName
    21         );
     13             "Usage: %s [-h] [-v] [-s NX NY] IN OUT\n"
     14             "where\n"
     15             "\t-h           Help (this info)\n"
     16             "\t-v           Increase verbosity level\n"
     17             "\t-s NX NY     Size of output image\n"
     18             "\tIN           Input image, which has associated .map file\n"
     19             "\tOUT          Output image\n",
     20             programName
     21        );
    2222}
    2323
    … …  
    3838
    3939    // Parameters with default values
    40     int outnx = 0, outny = 0;           // Size of output image
    41     int verbose = 0;                    // Verbosity level
     40    int outnx = 0, outny = 0;           // Size of output image
     41    int verbose = 0;                    // Verbosity level
     42    float bad = 0.0;                    // Bad level
    4243
    4344    // Parse command line
    44     const char *programName = argv[0];  // Program name
     45    const char *programName = argv[0];  // Program name
    4546    /* Variables for getopt */
    4647    int opt;   /* Option, from getopt */
    … …  
    4950
    5051    /* Parse command-line arguments using getopt */
    51     while ((opt = getopt(argc, argv, "hvs:")) != -1) {
     52    while ((opt = getopt(argc, argv, "hvb:s:")) != -1) {
    5253        switch (opt) {
    53           case 'h':
    54             help(programName);
    55             exit(EXIT_SUCCESS);
    56           case 'v':
     54          case 'h':
     55            help(programName);
     56            exit(EXIT_SUCCESS);
     57          case 'v':
    5758            verbose++;
    5859            break;
    59           case 's':
     60          case 's':
    6061            if ((sscanf(argv[optind-1], "%d", &outnx) != 1) || (sscanf(argv[optind++], "%d", &outny) != 1)) {
    6162                // Note: incrementing optind, so I can read more than one parameter.
    … …  
    6364                exit(EXIT_FAILURE);
    6465            }
    65             break;
    66           default:
    67             help(programName);
    68             exit(EXIT_FAILURE);
    69         }
     66            break;
     67          case 'b':
     68            if (sscanf(optarg, "%f", &bad) != 1) {
     69                help(programName);
     70                exit(EXIT_FAILURE);
     71            }
     72            break;
     73          default:
     74            help(programName);
     75            exit(EXIT_FAILURE);
     76        }
    7077    }
    7178
    … …  
    7784        exit(EXIT_FAILURE);
    7885    }
    79     const char *inName = argv[0];       // Input filename
    80     const char *outName = argv[1];      // Output filename
     86    const char *inName = argv[0];       // Input filename
     87    const char *outName = argv[1];      // Output filename
    8188
    82     psFits *imageFile = psFitsAlloc(inName);
    83     psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
    84     psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
     89    psFits *imageFile = psFitsOpen(inName, "r");
     90    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
     91    psMetadata *header = psFitsReadHeader(NULL, imageFile); // FITS header
     92    psImage *image = psFitsReadImage(NULL, imageFile, imageRegion, 0);
    8593    if (image == NULL) {
    86         psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", inName);
    87         exit(EXIT_FAILURE);
     94        psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", inName);
     95        exit(EXIT_FAILURE);
    8896    }
    8997    psTrace("shift", 4, "Image %s is %dx%d\n", inName, image->numCols, image->numRows);
    9098    // Convert to 32-bit floating point, in necessary
    9199    if (image->type.type != PS_TYPE_F32) {
    92         psTrace("shift", 5, "Converting %s to floating point in memory....", inName);
    93         psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
    94         psFree(image);
    95         image = temp;
     100        psTrace("shift", 5, "Converting %s to floating point in memory....", inName);
     101        psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
     102        psFree(image);
     103        image = temp;
    96104    }
    97     psFree(imageFile);
    98     psFree(imageRegion);
     105    psFitsClose(imageFile);
     106
     107    // Generate masks
     108    psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
     109    for (int y = 0; y < image->numRows; y++) {
     110        for (int x = 0; x < image->numCols; x++) {
     111            if (image->data.F32[y][x] <= bad) {
     112                mask->data.U8[y][x] = 1;
     113            } else {
     114                mask->data.U8[y][x] = 0;
     115            }
     116        }
     117    }
    99118
    100119    // Read map
    101     char mapName[MAXCHAR];              // Name of map file
     120    char mapName[MAXCHAR];              // Name of map file
    102121    sprintf(mapName, "%s.map", inName);
    103122    psPlaneTransform *map = stacReadMap(mapName);
    104123
    105124    // Functions work on array of images, so:
    106     psArray *images = psArrayAlloc(1); // An array of one
    107     psArray *maps = psArrayAlloc(1);// An array of one
     125    psArray *images = psArrayAlloc(1);  // An array of one
     126    psArray *maps = psArrayAlloc(1);    // An array of one
     127    psArray *masks = psArrayAlloc(1);   // An array of one
    108128    images->data[0] = image;
    109129    maps->data[0] = map;
     130    masks->data[0] = mask;
    110131
    111132    // Get size
    112133    if (outnx == 0 || outny == 0) {
    113         stacSize(&outnx, &outny, NULL, NULL, images, maps);
     134        stacSize(&outnx, &outny, NULL, NULL, images, maps);
    114135    }
    115136
    … …  
    119140    // Transform inputs and errors
    120141    psArray *transformed = NULL;
    121     (void)stacTransform(&transformed, NULL, images, inverseMaps, NULL, NULL, NULL, NULL, NULL, outnx, outny);
     142    (void)stacTransform(&transformed, NULL, images, inverseMaps, NULL, masks, NULL, NULL, NULL, outnx, outny);
     143
     144    // Update FITS header appropriately
     145    psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of pixels in x",
     146                     (int)((psImage*)transformed->data[0])->numCols);
     147    psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of pixels in y",
     148                     (int)((psImage*)transformed->data[0])->numRows);
     149    psMetadataAddS32(header, PS_LIST_TAIL, "BITPIX", PS_META_REPLACE, "Bits per pixel", -32);
    122150
    123151    // Write out transformed image
    124     psFits *outFile = psFitsAlloc(outName);
    125     if (!psFitsWriteImage(outFile, NULL, transformed->data[0], 0, NULL)) {
    126         psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
     152    psFits *outFile = psFitsOpen(outName, "w");
     153    int numPix = psImageClipNaN(transformed->data[0], 0.0);
     154    if (numPix > 0) {
     155        psTrace("stac", 3, "Clipping %d NaN pixels to zero.\n", numPix);
     156    }
     157    if (!psFitsWriteImage(outFile, header, transformed->data[0], 0)) {
     158        psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
    127159    }
    128160    psTrace("shift", 1, "Transformed image written to %s\n", outName);
    129     psFree(outFile);
     161    psFitsClose(outFile);
    130162
    131163    // Free everything I've used
    132164    psFree(images);
    133165    psFree(maps);
     166    psFree(masks);
    134167    psFree(inverseMaps);
    135168    psFree(transformed);

  • trunk/stac/src/shiftSize.c

    r3667 r5743  
    7979
    8080    // Read input files
    81     psArray *images = stacReadImages(inputs);
     81    psArray *images = stacReadImages(NULL, inputs);
    8282
    8383    // Read maps

  • trunk/stac/src/stac.c

    r3673 r5743  
    1818    // Set trace levels
    1919    psTraceSetLevel(".",0);
     20    psTraceSetLevel("stac",10);
    2021    psTraceSetLevel("stac.checkMemory",10);
    2122    psTraceSetLevel("stac.config",10);
    … …  
    3839
    3940    // Read input files
    40     psArray *inputs = stacReadImages(config->inputs);
     41    psArray *headers = NULL;            // Array of headers
     42    psArray *inputs = stacReadImages(&headers, config->inputs);
     43
     44    // Generate masks
     45    psArray *masks = psArrayAlloc(inputs->n);
     46    for (int i = 0; i < inputs->n; i++) {
     47        psImage *image = inputs->data[i]; // Image for which to get mask
     48        psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
     49        for (int y = 0; y < image->numRows; y++) {
     50            for (int x = 0; x < image->numCols; x++) {
     51                if (image->data.F32[y][x] <= config->bad) {
     52                    mask->data.U8[y][x] = 1;
     53                } else {
     54                    mask->data.U8[y][x] = 0;
     55                }
     56            }
     57        }
     58        masks->data[i] = mask;
     59    }
    4160
    4261    // Read maps
    … …  
    4766    // Get size, if not input
    4867    if (config->outnx == 0 || config->outny == 0) {
    49         stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps);
     68        stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps);
     69    }
     70
     71    // Fix up the headers
     72    for (int i = 0; i < headers->n; i++) {
     73        psMetadata *header = headers->data[i];
     74        psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of pixels in x",
     75                         config->outnx);
     76        psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of pixels in y",
     77                         (int)config->outny);
     78        psMetadataAddS32(header, PS_LIST_TAIL, "BITPIX", PS_META_REPLACE, "Bits per pixel", -32);
     79#if 0
     80        bool mdok = true;               // Result of MD lookup
     81        float crpix1 = psMetadataLookupF32(&mdok, header, "CRPIX1");
     82        if (mdok && !isnan(crpix1)) {
     83            psMetadataAddF32(header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE,
     84                             "WCS Coordinate reference pixel", crpix1 -
     85#endif
    5086    }
    5187
    … …  
    5894    // Write error images out to check
    5995    for (int i = 0; i < inputs->n; i++) {
    60         char errName[MAXCHAR];          // Filename of error image
    61         sprintf(errName,"%s.err",config->inputs->data[i]);
    62         psFits *errorFile = psFitsAlloc(errName);
    63         if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0, NULL)) {
    64             psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
    65         }
    66         psTrace("stac", 1, "Error image written to %s\n", errName);
    67         psFree(errorFile);
     96        char errName[MAXCHAR];          // Filename of error image
     97        sprintf(errName,"%s.err",config->inputs->data[i]);
     98        psFits *errorFile = psFitsOpen(errName, "w");
     99        if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0)) {
     100            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
     101        }
     102        psTrace("stac", 1, "Error image written to %s\n", errName);
     103        psFitsClose(errorFile);
    68104    }
    69105#endif
    … …  
    72108    psArray *transformed = NULL;
    73109    psArray *transformedErrors = NULL;
    74     (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL,
    75                         config->outnx, config->outny);
     110    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks, NULL, NULL, NULL,
     111                        config->outnx, config->outny);
    76112    psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime);
    77113
    … …  
    79115    // Write transformed images out to check
    80116    for (int i = 0; i < inputs->n; i++) {
    81         char shiftName[MAXCHAR];        // Filename of shift image
    82         char errName[MAXCHAR];          // Filename of error image
    83         sprintf(shiftName,"%s.shift1",config->inputs->data[i]);
    84         sprintf(errName,"%s.shifterr1",config->inputs->data[i]);
    85         psFits *shiftFile = psFitsAlloc(shiftName);
    86         psFits *errFile = psFitsAlloc(errName);
    87         psImage *trans = transformed->data[i]; // Transformed image
    88         psImage *transErr = transformedErrors->data[i]; // Transformed error image
    89         if (!psFitsWriteImage(shiftFile, NULL, trans, 0, NULL)) {
    90             psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
    91         }
    92         psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
    93         if (!psFitsWriteImage(errFile, NULL, transErr, 0, NULL)) {
    94             psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
    95         }
    96         psTrace("stac", 1, "Shifted error image written to %s\n", errName);
    97         psFree(shiftFile);
    98         psFree(errFile);
    99     }
    100 #endif
    101 
    102 
    103 
    104     psVector *scales = NULL;            // Relative scales between images
    105     psVector *offsets = NULL;           // Offsets between images
     117        char shiftName[MAXCHAR];        // Filename of shift image
     118        char errName[MAXCHAR];          // Filename of error image
     119        sprintf(shiftName,"%s.shift1",config->inputs->data[i]);
     120        sprintf(errName,"%s.shifterr1",config->inputs->data[i]);
     121        psFits *shiftFile = psFitsOpen(shiftName, "w");
     122        psFits *errFile = psFitsOpen(errName, "w");
     123        psImage *trans = transformed->data[i]; // Transformed image
     124        psImage *transErr = transformedErrors->data[i]; // Transformed error image
     125        if (!psFitsWriteImage(shiftFile, NULL, trans, 0)) {
     126            psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
     127        }
     128        psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
     129        if (!psFitsWriteImage(errFile, NULL, transErr, 0)) {
     130            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
     131        }
     132        psTrace("stac", 1, "Shifted error image written to %s\n", errName);
     133        psFitsClose(shiftFile);
     134        psFitsClose(errFile);
     135    }
     136#endif
     137
     138    psVector *scales = NULL;            // Relative scales between images
     139    psVector *offsets = NULL;           // Offsets between images
    106140    (void)stacScales(&scales, &offsets, transformed, config->starFile, config->starMapFile, config->xMapDiff,
    107                      config->yMapDiff, config->aper);
     141                     config->yMapDiff, config->aper);
    108142    // Rescale the images
    109143    (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets);
    110144    // Set the saturation and bad values
    111145    psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits
    112     psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits
     146    psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits
    113147    for (int i = 0; i < transformed->n; i++) {
    114         saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
    115         bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
     148        saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
     149        bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
     150    }
     151
     152    // Save shifted images
     153    if (config->saveShifts) {
     154        psImage *image = NULL;          // Copy of image, to remove NAN for output
     155        for (int i = 0; i < transformed->n; i++) {
     156            char shiftName[MAXCHAR];    // Filename of shifted image
     157            sprintf(shiftName, "%s.shift", config->inputs->data[i]);
     158            psFits *shiftFile = psFitsOpen(shiftName, "w");
     159            image = psImageCopy(NULL, transformed->data[i], PS_TYPE_F32);
     160            (void)psImageClipNaN(image, 0.0);
     161            if (!psFitsWriteImage(shiftFile, headers->data[i], image, 0)) {
     162                psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
     163            }
     164            psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName);
     165            psFitsClose(shiftFile);
     166            psFree(image);
     167        }
    116168    }
    117169
    … …  
    120172    psImage *combined = NULL;
    121173    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, saturated,
    122                       bad, config->reject);
     174                      bad, config->reject);
    123175
    124176    psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime);
    … …  
    128180    // Write rejection images out to check
    129181    for (int i = 0; i < rejected->n; i++) {
    130         char rejName[MAXCHAR];  // Filename of rejection image
    131         sprintf(rejName, "%s.shiftrej", config->inputs->data[i]);
    132        
    133         psFits *rejFile = psFitsAlloc(rejName);
    134         if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0, NULL)) {
    135             psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
    136         }
    137         psTrace("stac", 1, "Rejection image written to %s\n", rejName);
    138         psFree(rejFile);
     182        char rejName[MAXCHAR];  // Filename of rejection image
     183        sprintf(rejName, "%s.shiftrej", config->inputs->data[i]);
     184
     185        psFits *rejFile = psFitsOpen(rejName, "w");
     186        if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0)) {
     187            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
     188        }
     189        psTrace("stac", 1, "Rejection image written to %s\n", rejName);
     190        psFitsClose(rejFile);
    139191    }
    140192
    141193    // Write out pre-combined image
    142     char preName[MAXCHAR];              // Filename of precombined image
     194    char preName[MAXCHAR];              // Filename of precombined image
    143195    sprintf(preName, "%s.pre", config->output);
    144     psFits *preFile = psFitsAlloc(preName);
    145     if (!psFitsWriteImage(preFile, NULL, combined, 0, NULL)) {
    146         psErrorStackPrint(stderr, "Unable to write image: %s\n", preName);
     196    psFits *preFile = psFitsOpen(preName, "w");
     197    if (!psFitsWriteImage(preFile, NULL, combined, 0)) {
     198        psErrorStackPrint(stderr, "Unable to write image: %s\n", preName);
    147199    }
    148200    psTrace("stac", 1, "Pre-combined image written to %s\n", preName);
    149     psFree(preFile);
     201    psFitsClose(preFile);
    150202#endif
    151203
    … …  
    153205    psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest
    154206    for (int i = 0; i < inputs->n; i++) {
    155         regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
    156                                               ((psImage*)(inputs->data[i]))->numCols,
    157                                               ((psImage*)(inputs->data[i]))->numRows);
    158 #ifdef TESTING
    159         char regionName[MAXCHAR];       // Filename of region image
    160         sprintf(regionName,"%s.region",config->inputs->data[i]);
    161         psFits *regionFile = psFitsAlloc(regionName);
    162         if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0, NULL)) {
    163             psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName);
    164         }
    165         psTrace("stac", 1, "Region image written to %s\n", regionName);
    166         psFree(regionFile);
     207        regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
     208                                              ((psImage*)(inputs->data[i]))->numCols,
     209                                              ((psImage*)(inputs->data[i]))->numRows);
     210#ifdef TESTING
     211        char regionName[MAXCHAR];       // Filename of region image
     212        sprintf(regionName,"%s.region",config->inputs->data[i]);
     213        psFits *regionFile = psFitsOpen(regionName, "w");
     214        if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0)) {
     215            psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName);
     216        }
     217        psTrace("stac", 1, "Region image written to %s\n", regionName);
     218        psFitsClose(regionFile);
    167219#endif
    168220    }
    … …  
    170222    // Transform rejected pixels to source frame
    171223    psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config->frac,
    172                                             config->grad, config->inputs);
     224                                            config->grad, config->inputs);
    173225
    174226    // Get regions of interest in the output frame
    175227    psImage *combineRegion = psImageAlloc(config->outnx, config->outny, PS_TYPE_U8);
    176228    for (int y = 0; y < config->outny; y++) {
    177         for (int x = 0; x < config->outnx; x++) {
    178             combineRegion->data.U8[y][x] = 0;
    179         }
    180     }
    181     for (int i = 0; i < inputs->n; i++) {
    182         psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx,
    183                                              config->outny);
    184         psBinaryOp(combineRegion, combineRegion, "+", region);
    185         psFree(region);
     229        for (int x = 0; x < config->outnx; x++) {
     230            combineRegion->data.U8[y][x] = 0;
     231        }
     232    }
     233    for (int i = 0; i < inputs->n; i++) {
     234        psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx,
     235                                             config->outny);
     236        psBinaryOp(combineRegion, combineRegion, "+", region);
     237        psFree(region);
     238
     239        // Do OR of masks
     240        psImage *mask = masks->data[i]; // Mask of input image (bad columns etc)
     241        psImage *cr = rejectedSource->data[i]; // Mask of CRs
     242        for (int y = 0; y < mask->numRows; y++) {
     243            for (int x = 0; x < mask->numCols; x++) {
     244                if (cr->data.U8[y][x] > 0) {
     245                    mask->data.U8[y][x] = 1;
     246                }
     247            }
     248        }
    186249    }
    187250
    188251    // Redo transformation with the masks and scales/offsets
    189     (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,
    190                         combineRegion, scales, offsets, config->outnx, config->outny);
     252    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks,
     253                        combineRegion, scales, offsets, config->outnx, config->outny);
    191254
    192255    // Combine the newly-transformed CR-free images, no rejection
    … …  
    194257    rejected = NULL;
    195258    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, saturated,
    196                       bad, config->reject);
     259                      bad, config->reject);
    197260
    198261    // Write output image
    199     psFits *outFile = psFitsAlloc(config->output);
    200     if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) {
    201         psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output);
     262    psFits *outFile = psFitsOpen(config->output, "w");
     263    int numPix = psImageClipNaN(combined, 0.0);
     264    if (numPix > 0) {
     265        psTrace("stac", 3, "Clipping %d NaN pixels to zero.\n", numPix);
     266    }
     267    if (!psFitsWriteImage(outFile, headers->data[0], combined, 0)) {
     268        psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output);
    202269    }
    203270    psTrace("stac", 1, "Combined image written to %s\n", config->output);
    204     psFree(outFile);
    205 
    206     // Save shifted images
    207     if (config->saveShifts) {
    208         for (int i = 0; i < transformed->n; i++) {
    209             char shiftName[MAXCHAR];    // Filename of shifted image
    210             sprintf(shiftName, "%s.shift", config->inputs->data[i]);
    211             psFits *shiftFile = psFitsAlloc(shiftName);
    212             if (!psFitsWriteImage(shiftFile, NULL, transformed->data[i], 0, NULL)) {
    213                 psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
    214             }
    215             psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName);
    216             psFree(shiftFile);
    217         }
    218     }
     271    psFitsClose(outFile);
    219272
    220273    // Free everything I've used
    … …  
    228281    psFree(inverseMaps);
    229282    psFree(errors);
     283    psFree(masks);
    230284    psFree(transformedErrors);
    231285    psFree(transformed);
    232286    psFree(rejectedSource);
    233287    psFree(combined);
     288    psFree(saturated);
     289    psFree(bad);
    234290
    235291    psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime);

  • trunk/stac/src/stac.h

    r3673 r5743  
    1818
    1919// Read the input files and return an array of images
    20 psArray *stacReadImages(psArray *filenames // The file names of the images
     20psArray *stacReadImages(psArray **headers, // The image headers, to be returned
     21                        psArray *filenames // The file names of the images
    2122    );
    2223
    … …  
    7273
    7374// Print out a memblock when it's allocated --- this function used as a callback
    74 psMemoryId stacMemPrint(const psMemBlock *ptr);
     75psMemId stacMemPrint(const psMemBlock *ptr);
    7576
    7677

  • trunk/stac/src/stacAreaOfInterest.c

    r3375 r5743  
    1313    )
    1414{
    15     psDPolynomial2D *xPoly = transform->x;
    16     psDPolynomial2D *yPoly = transform->y;
     15    psPolynomial2D *xPoly = transform->x;
     16    psPolynomial2D *yPoly = transform->y;
    1717
    1818    assert(xPoly->type == PS_POLYNOMIAL_ORD);

  • trunk/stac/src/stacCheckMemory.c

    r3375 r5743  
    77
    88
    9 psMemoryId stacMemPrint(const psMemBlock *ptr)
     9psMemId stacMemPrint(const psMemBlock *ptr)
    1010{
    1111    psLogMsg("stac.memoryPrint", PS_LOG_INFO,
    … …  
    4141    psTrace("stac.checkMemory", 1, "Checking for memory problems....\n");
    4242
    43     (void)psMemProblemCallbackSet(stacMemoryProblem); // Set callback for corruption
     43    (void)psMemProblemCallbackSet((psMemProblemCallback)stacMemoryProblem); // Set callback for corruption
    4444
    4545    if ((leakFile = fopen(LEAKS, "w")) == NULL) {

  • trunk/stac/src/stacCombine.c

    r3673 r5743  
    1010float stacCombineMean(psVector *values, // Values for which to take the mean
    1111                      psVector *errors, // Errors in the values
    12                       psVector *masks   // Masks for the values, 0 = don't use, 1 = use
     12                      psVector *masks   // Masks for the values, 1 = don't use, 0 = use
    1313    )
    1414{
    … …  
    149149                    }
    150150                }
    151            
     151               
    152152                float average = stacCombineMean(pixels, deltas, mask); // Combined value
     153
     154                // We set the value BEFORE the rejection iteration because not all pixels that we reject
     155                // here will be rejected in the final cut.
     156                (*combined)->data.F32[y][x] = average;
    153157
    154158#ifdef TESTING
    … …  
    188192                }
    189193           
    190                 (*combined)->data.F32[y][x] = average;
    191 
    192194            } // Pixels of interest
    193195
    … …  
    201203    sprintf(chiName,"chi2_%d.fits",iteration);
    202204    psFits *chiFile = psFitsAlloc(chiName);
    203     if (!psFitsWriteImage(chiFile, NULL, chi2 , 0, NULL)) {
     205    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0)) {
    204206        psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
    205207    }

  • trunk/stac/src/stacConfig.c

    r3666 r5743  
    2626    config->saturated = 65535.0;        // Saturation level
    2727    config->bad = 0.0;                  // Bad level
    28     config->reject = 3.0;
    29     config->frac = 0.5;
    30     config->grad = 0.6;
     28    config->reject = 2.5;
     29    config->frac = 0.45;
     30    config->grad = 0.7;
    3131    config->nReject = 2;
    3232

  • trunk/stac/src/stacInvertMaps.c

    r3666 r5743  
    2727        assert(oldMap->x->nX == oldMap->x->nY && oldMap->y->nX == oldMap->y->nY &&
    2828               oldMap->x->nX == oldMap->y->nX);
    29         int order = oldMap->x->nX;      // Polynomial order
     29        int order = oldMap->x->nX + 1;  // Polynomial order
    3030        psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
    31         psPlaneTransform *newMap = psPlaneTransformAlloc(order, order); // Inverted map
     31        psPlaneTransform *newMap = psPlaneTransformAlloc(order + 1, order + 1); // Inverted map
    3232
    3333        // Create fake polynomial to use in evaluation
    34         psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
     34        psPolynomial2D *fakePoly = psPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
    3535        for (int i = 0; i < order; i++) {
    3636            for (int j = 0; j < order; j++) {
    … …  
    8282
    8383                    fakePoly->mask[i][j] = 0;
    84                     double ijPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
     84                    double ijPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
    8585                    fakePoly->mask[i][j] = 1;
    8686
    … …  
    8989                           
    9090                            fakePoly->mask[m][n] = 0;
    91                             double mnPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
     91                            double mnPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
    9292                            fakePoly->mask[m][n] = 1;
    9393

  • trunk/stac/src/stacRead.c

    r3681 r5743  
    11#include <stdio.h>
    22#include <string.h>
     3#include <assert.h>
    34#include "pslib.h"
    45#include "stac.h"
    56
    6 #define BUFFER 100                      // Size of buffer for incrementally reading coordinates
    7 
    8 psArray *stacReadImages(psArray *filenames // The file names of the images
     7#define BUFFER 100                      // Size of buffer for incrementally reading coordinates
     8
     9psArray *stacReadImages(psArray **headers, // The image headers, to be returned
     10    psArray *filenames // The file names of the images
    911    )
    1012{
    11     int nFiles = filenames->n;          // The number of input files
     13    int nFiles = filenames->n;          // The number of input files
    1214    psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
    13     psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
     15    assert(!headers || ! *headers || (*headers)->n == nFiles);
     16    if (headers && ! *headers) {
     17        *headers = psArrayAlloc(nFiles);
     18    }
     19
     20    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
    1421
    1522    psTrace("stac.read.images", 1, "Reading input images....\n");
    1623    for (int i = 0; i < nFiles; i++) {
    17         psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
    18         psFits *imageFile = psFitsAlloc(filenames->data[i]);
    19         // We only read PHUs --- not mucking around with extensions for now
    20         psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
    21         if (image == NULL) {
    22             psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
    23             exit(EXIT_FAILURE);
    24         }
    25         psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
    26                 image->numRows);
    27         // Convert to 32-bit floating point, in necessary
    28         if (image->type.type != PS_TYPE_F32) {
    29             psTrace("stac.read.images", 3, "Converting %s to floating point in memory....",
    30                     filenames->data[i]);
    31             images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
    32             psFree(image);
    33         } else {
    34             int numNaN = psImageClipNaN(image, FLT_MAX);
    35             if (numNaN) {
    36                 psTrace("stac.read.images", 5, "Clipped %d NaN pixels.\n", numNaN);
    37             }
    38             images->data[i] = image;
    39         }
    40         psFree(imageFile);
     24        psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
     25        psFits *imageFile = psFitsOpen(filenames->data[i], "r");
     26        // We only read PHUs --- not mucking around with extensions for now
     27        if (headers) {
     28            (*headers)->data[i] = psFitsReadHeader(NULL, imageFile);
     29        }
     30        psImage *image = psFitsReadImage(NULL, imageFile, imageRegion, 0);
     31        if (image == NULL) {
     32            psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
     33            exit(EXIT_FAILURE);
     34        }
     35        psFitsClose(imageFile);
     36        psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
     37                image->numRows);
     38        // Convert to 32-bit floating point, in necessary
     39        if (image->type.type != PS_TYPE_F32) {
     40            psTrace("stac.read.images", 3, "Converting %s to floating point in memory....\n",
     41                    filenames->data[i]);
     42            images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
     43            psFree(image);
     44        }
     45        int numNaN = psImageClipNaN(image, 0.0);
     46        if (numNaN) {
     47            psTrace("stac.read.images", 5, "Clipped %d NaN pixels to zero.\n", numNaN);
     48        }
     49        images->data[i] = image;
    4150    }
    4251    psTrace("stac.read.images",1,"%d input images read.\n",nFiles);
    43     psFree(imageRegion);
    4452
    4553    return images;
    … …  
    5058    FILE *file = fopen(filename, "r");
    5159    if (file == NULL) {
    52         psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
    53         return NULL;
     60        psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
     61        return NULL;
    5462    }
    5563
    … …  
    5765
    5866    psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
    59     float x, y;                         // Coordinates to read
    60     int num = 0;                        // Number of coordinates read
     67    float x, y;                         // Coordinates to read
     68    int num = 0;                        // Number of coordinates read
    6169    while (fscanf(file, "%f %f\n", &x, &y) == 2) {
    62         psPlane *coord = psPlaneAlloc();// A coordinate
    63         coord->x = x;
    64         coord->y = y;
    65         coords->data[num] = coord;
    66         num++;
    67         if (num % BUFFER) {
    68             coords = psArrayRealloc(coords, num + BUFFER);
    69         }
     70        psPlane *coord = psPlaneAlloc();// A coordinate
     71        coord->x = x;
     72        coord->y = y;
     73        coords->data[num] = coord;
     74        num++;
     75        if (num % BUFFER) {
     76            coords = psArrayRealloc(coords, num + BUFFER);
     77        }
    7078    }
    7179    coords->n = num;
    … …  
    8290    FILE *mapfp = fopen(filename, "r");
    8391    if (mapfp == NULL) {
    84         psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
    85         return NULL;
     92        psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
     93        return NULL;
    8694    }
    8795    // Read the file
    8896    psTrace("stac.read.map", 5, "Reading map file %s....\n", filename);
    89    
     97
    9098    // Format is now:
    9199    // order
    92100    // x coefficients
    93101    // y coefficients
    94     // 
     102    //
    95103    // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
    96104    // and the coefficients are read by the following pseudo-code:
    97     // 
    98     //  for (int k = 0; k < order + 1; k++)
    99     //      for (int j = 0; j < k; j++)
    100     //          int i = k - j;
     105    //
     106    //  for (int k = 0; k < order + 1; k++)
     107    //      for (int j = 0; j < k; j++)
     108    //          int i = k - j;
    101109    //              read coefficient of x^i y^j
    102     // 
     110    //
    103111    // This produces the ordering:
    104112    // x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
    105     // 
     113    //
    106114    // This is, of course, for third order polynomials.
    107115    // For lower orders, the list is truncated at the appropriate level.
    108116
    109     int order = 0;                      // Polynomial order
     117    int order = 0;                      // Polynomial order
    110118    if (fscanf(mapfp, "%d", &order) != 1) {
    111         psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
    112         fclose(mapfp);
    113         return NULL;
    114     }
    115    
     119        psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     120        fclose(mapfp);
     121        return NULL;
     122    }
     123
    116124    psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);
    117    
     125
    118126    psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
    119127    // Set coefficient masks
    120128    for (int i = 0; i < order + 1; i++) {
    121         for (int j = 0; j < order + 1; j++) {
    122             if (i + j > order + 1) {
    123                 map->x->mask[i][j] = 1;
    124                 map->y->mask[i][j] = 1;
    125             } else {
    126                 map->x->mask[i][j] = 0;
    127                 map->y->mask[i][j] = 0;
    128             }
    129         }
    130     }
    131    
     129        for (int j = 0; j < order + 1; j++) {
     130            if (i + j > order + 1) {
     131                map->x->mask[i][j] = 1;
     132                map->y->mask[i][j] = 1;
     133            } else {
     134                map->x->mask[i][j] = 0;
     135                map->y->mask[i][j] = 0;
     136            }
     137        }
     138    }
     139
    132140    // Read x coefficients
    133141    psTrace("stac.read.map", 7, "x' = \n");
    134142    for (int k = 0; k < order + 1; k++) {
    135         for (int j = 0; j < k + 1; j++) {
    136             int i = k - j;
    137             if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
    138                 psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
    139                 fclose(mapfp);
    140                 psFree(map);
    141                 return NULL;
    142             }
    143             psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
    144         }
     143        for (int j = 0; j < k + 1; j++) {
     144            int i = k - j;
     145            if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
     146                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     147                fclose(mapfp);
     148                psFree(map);
     149                return NULL;
     150            }
     151            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
     152        }
    145153    }
    146154    // Read y coefficients
    147155    psTrace("stac.read.maps", 7, "y' = \n");
    148156    for (int k = 0; k < order + 1; k++) {
    149         for (int j = 0; j < k + 1; j++) {
    150             int i = k - j;
    151             if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
    152                 psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
    153                 fclose(mapfp);
    154                 psFree(map);
    155                 return NULL;
    156             }
    157             psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
    158         }
    159     }
    160    
     157        for (int j = 0; j < k + 1; j++) {
     158            int i = k - j;
     159            if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
     160                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     161                fclose(mapfp);
     162                psFree(map);
     163                return NULL;
     164            }
     165            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
     166        }
     167    }
     168
    161169    fclose(mapfp);
    162170
    163171    return map;
    164172}
    165    
     173
    166174
    167175
    … …  
    169177    )
    170178{
    171     int nFiles = filenames->n;          // The number of input files
     179    int nFiles = filenames->n;          // The number of input files
    172180    psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
    173     char mapfile[MAXCHAR];              // Filename of map
    174    
     181    char mapfile[MAXCHAR];              // Filename of map
     182
    175183    psTrace("stac.read.maps",1,"Reading maps....\n");
    176184    for (int i = 0; i < nFiles; i++) {
    177         if (strlen(filenames->data[i]) > MAXCHAR - 4) {
    178             psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
    179             exit(EXIT_FAILURE);
    180         }
    181         // Read the file
    182         sprintf(mapfile,"%s.map",filenames->data[i]);
    183         maps->data[i] = stacReadMap(mapfile);
    184         if (maps->data[i] == NULL) {
    185             psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
    186         }
     185        if (strlen(filenames->data[i]) > MAXCHAR - 4) {
     186            psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
     187            exit(EXIT_FAILURE);
     188        }
     189        // Read the file
     190        sprintf(mapfile,"%s.map",filenames->data[i]);
     191        maps->data[i] = stacReadMap(mapfile);
     192        if (maps->data[i] == NULL) {
     193            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
     194        }
    187195
    188196    }

  • trunk/stac/src/stacRejection.c

    r3667 r5743  
    7474    psTrace("stac.rejection", 1, "Mapping rejection masks back to source....\n");
    7575
     76    // Vectors for calculating the mean gradient
     77    psVector *grads = psVectorAlloc(nImages, PS_TYPE_F32); // Gradient for each image
     78    psVector *gradsMask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask for gradient vector
     79
    7680    // Transform rejection masks back to source
    7781    psArray *inputRej = psArrayAlloc(nImages);
    … …  
    107111        // calculate derivatives of the map, and use that as a buffer around the transformed position
    108112        // in the input image.
     113        psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    109114        for (int y = 0; y < nyInput; y++) {
    110115            for (int x = 0; x < nxInput; x++) {
    … …  
    135140                                    (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
    136141                                    // Calculate the gradient
    137                                     meanGrads += stacGradient(inputs->data[j], xPix, yPix);
     142                                    grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
     143                                    gradsMask->data.U8[j] = 0;
    138144                                    numGrads++;
     145                                } else {
     146                                    gradsMask->data.U8[j] = 1; // Mask this one
    139147                                }
     148                            } else {
     149                                gradsMask->data.U8[j] = 1; // Mask this one
    140150                            }
    141151                        }
    142152                        if (numGrads > 0) {
    143                             meanGrads /= (float)numGrads;
     153                            (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
     154                            meanGrads = median->sampleMedian;
    144155                        } else {
    145                             meanGrads = 0;
     156                            meanGrads = 0.0;
    146157                        }
    147158
    148159#ifdef TESTING
    149                         gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
    150 #endif
    151 
    152                         if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
     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) {
    153166                            mask->data.U8[y][x] = 1;
    154167                            nBad++;
    … …  
    173186            }
    174187        } // Iterating over pixels
     188        psFree(median);
    175189
    176190#ifdef CRFLUX
    … …  
    194208        psFits *rejmapFile = psFitsAlloc(rejmapName);
    195209        psFits *gradFile = psFitsAlloc(gradName);
    196         if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
     210        if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
    197211            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
    198212        }
    199213        psTrace("stac", 1, "Mask image written to %s\n", maskName);
    200         if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0, NULL)) {
     214        if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0)) {
    201215            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
    202216        }
    203217        psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
    204         if (!psFitsWriteImage(gradFile, NULL, gradient, 0, NULL)) {
     218        if (!psFitsWriteImage(gradFile, NULL, gradient, 0)) {
    205219            psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
    206220        }
    … …  
    218232    }
    219233
     234    psFree(grads);
     235    psFree(gradsMask);
    220236
    221237    psFree(inCoords);

  • trunk/stac/src/stacScales.c

    r3673 r5743  
    1414    assert(image->type.type == PS_TYPE_F32);
    1515
    16 // Will use robust median instead of sampling --- it's supposed to be fast.
    1716#if 1
    1817    int size = image->numCols * image->numRows; // Number of pixels in image
    1918    int numSamples = size / sample;     // Number of samples in image
    2019    psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample
     20    psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
    2121
    2222    int offset = 0;                     // Offset from start of the row
    … …  
    2727        while (col < image->numCols) {
    2828            values->data.F32[index] = image->data.F32[row][col];
     29            if (isnan(values->data.F32[index])) {
     30                mask->data.U8[index] = 1;
     31            } else {
     32                mask->data.U8[index] = 0;
     33            }
    2934            col += sample;
    3035            index++;
    … …  
    3439
    3540    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    36     stats = psVectorStats(stats, values, NULL, NULL, 0);
     41    stats = psVectorStats(stats, values, NULL, mask, 1);
    3742    float median = stats->sampleMedian;
    3843    psFree(stats);
    3944    psFree(values);
     45    psFree(mask);
    4046#else
     47    // Will use robust median instead of sampling --- it's supposed to be fast.
    4148    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Using a clipped mean because median is SLOW
    4249    stats->clipSigma = 3.0;
    … …  
    162169                        }
    163170                    }
    164                     psTrace("stac.scales", 9, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
    165171                    // Subtract background, renormalise to account for circular aperture
    166172                    if (numPix > 0 && sum > 0) {
    167173                        sum -= offsets->data.F32[i] * (float)numPix;
    168174                        photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix;
     175                        if (photometry->data.F32[j] > 0 && finite(photometry->data.F32[j])) {
     176                            mask->data.U8[j] = 1;
     177                            psTrace("stac.scales", 8, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
     178                        } else {
     179                            mask->data.U8[j] = 0;
     180                        }
     181                    } else {
    169182                        mask->data.U8[j] = 0;
    170                     } else {
    171                         mask->data.U8[j] = 1;
    172183                    }
    173184                }
    … …  
    182193        psVector *ref = stars->data[0]; // The reference photometry
    183194        psVector *refMask = masks->data[0]; // The reference mask
     195#if 0
    184196        psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics
    185197        stats->clipSigma = 2.5;
    186198        stats->clipIter = 3;
     199#else
     200        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
     201#endif
    187202        for (int i = 1; i < scales->n; i++) {
    188             psVector *compare = stars->data[i]; // The comparison photometry
    189             psVector *compareMask = masks->data[i]; // The comparison mask
    190 
    191             compare = (psVector*)psBinaryOp(compare, compare, "/", ref);
    192             compareMask = (psVector*)psBinaryOp(compareMask, compareMask, "+", refMask);
    193 
    194             stats = psVectorStats(stats, compare, NULL, compareMask, 3); // Use maskVal of 3 to catch 1 and 2
    195 
     203            psVector *input = stars->data[i];   // The comparison photometry
     204            psVector *inputMask = masks->data[i]; // The comparison mask
     205
     206            psVector *compare = (psVector*)psBinaryOp(NULL, input, "/", ref);
     207            psVector *compareMask = (psVector*)psBinaryOp(NULL, inputMask, "*", refMask);
     208            (void)psBinaryOp(compareMask, psScalarAlloc(1, PS_TYPE_U8), "-", compareMask);
     209
     210#if 0
     211            psTrace("stac.scales", 9, "Getting scale for image %d...\n", i);
     212            for (int j = 0; j < compare->n; j++) {
     213                if (compareMask->data.U8[j] & 1) {
     214                    psTrace("stac.scales", 9, "Bad star %d: %f %f --> %f %d\n", j, input->data.F32[j],
     215                            ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
     216                } else {
     217                    psTrace("stac.scales", 9, "Good star %d: %f %f --> %f %d\n", j, input->data.F32[j],
     218                            ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
     219                }
     220            }
     221#endif
     222
     223            psVectorStats(stats, compare, NULL, compareMask, 1);
     224
     225#if 0
    196226            scales->data.F32[i] = stats->clippedMean;
     227#else
     228            scales->data.F32[i] = stats->sampleMedian;
     229#endif
    197230            psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
     231            psFree(compare);
     232            psFree(compareMask);
    198233        }
    199234        psFree(stats);

  • trunk/stac/src/stacTransform.c

    r3669 r5743  
    130130        for (int i = 0; i < nImages; i++) {
    131131            (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
     132            psImageInit((*outputs)->data[i], 0.0);
    132133        }
    133134    }
    … …  
    200201                    // Change PS_INTERPOLATE_BILINEAR to best available technique.
    201202                    outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
    202                                                                               detector->y, mask, 1, 0.0,
     203                                                                              detector->y, mask, 1, NAN,
    203204                                                                              PS_INTERPOLATE_BILINEAR);
    204205                    if (error) {
    … …  
    207208                                                                                                detector->x,
    208209                                                                                                detector->y,
    209                                                                                                 mask, 1, 0.0);
     210                                                                                                mask, 1, NAN);
    210211                    }
    211212

  • trunk/stac/src/stacWrite.c

    r3680 r5743  
    1616    }
    1717
    18     psDPolynomial2D *xMap = map->x;     // x transform
    19     psDPolynomial2D *yMap = map->y;     // y transform
     18    psPolynomial2D *xMap = map->x;      // x transform
     19    psPolynomial2D *yMap = map->y;      // y transform
    2020
    2121    // A crucial limitation of the current system --- the order of each polynomial must be the same
    2222    assert(xMap->nX == xMap->nY && yMap->nX == yMap->nY && xMap->nX == yMap->nX);
    23     int order = xMap->nX - 1;   // The polynomial order
     23    int order = xMap->nX;       // The polynomial order
    2424    fprintf(mapFile, "%d\n", order);
    2525   
Note: See TracChangeset for help on using the changeset viewer.