IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5742


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

Updating to use pslib-0.9.0; replaced psFitsAlloc with psFitsOpen

Location:
trunk/pois/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/pois/src/pois.c

    r5717 r5742  
    1818{
    1919    struct timeval tv;
    20  
     20
    2121    gettimeofday(&tv, NULL);
    2222    return(tv.tv_sec + 1.e-6 * tv.tv_usec);
    … …  
    3535
    3636    if (config->verbose) {
    37         psTraceSetLevel("pois", 10);
    38         psTraceSetLevel("pois.config", 10);
    39         psTraceSetLevel("pois.time", 10);
    40         psTraceSetLevel("pois.solution", 0);
    41         psTraceSetLevel("pois.kernelBasisFunctions", 10);
    42         psTraceSetLevel("pois.calculateEquation", 10);
    43         psTraceSetLevel("pois.convolveImage", 10);
    44         psTraceSetLevel("pois.findStamps", 10);
    45         psTraceSetLevel("pois.parseConfig", 10);
    46         psTraceSetLevel("pois.makeMask", 10);
    47         psTraceSetLevel("pois.extractKernel", 10);
    48         psTraceSetLevel("pois.calculateDeviations", 10);
    49         psTraceSetLevel("pois.checkKernel", 10);
    50         if (config->stampFile) {
    51             psTraceSetLevel("pois.checkStamp", 10);
    52         } else {
    53             psTraceSetLevel("pois.checkStamp", 0);
    54         }
    55         // Set logging level
    56         psLogSetLevel(9);
     37        psTraceSetLevel("pois", 10);
     38        psTraceSetLevel("pois.config", 10);
     39        psTraceSetLevel("pois.time", 10);
     40        psTraceSetLevel("pois.solution", 0);
     41        psTraceSetLevel("pois.kernelBasisFunctions", 10);
     42        psTraceSetLevel("pois.calculateEquation", 10);
     43        psTraceSetLevel("pois.convolveImage", 10);
     44        psTraceSetLevel("pois.findStamps", 10);
     45        psTraceSetLevel("pois.parseConfig", 10);
     46        psTraceSetLevel("pois.makeMask", 10);
     47        psTraceSetLevel("pois.extractKernel", 10);
     48        psTraceSetLevel("pois.calculateDeviations", 10);
     49        psTraceSetLevel("pois.checkKernel", 10);
     50        if (config->stampFile) {
     51            psTraceSetLevel("pois.checkStamp", 10);
     52        } else {
     53            psTraceSetLevel("pois.checkStamp", 0);
     54        }
     55        // Set logging level
     56        psLogSetLevel(9);
    5757    }
    5858
    5959    // Read inputs
    60     psMetadata *header = NULL;          // Header for output image
     60    psMetadata *header = NULL;          // Header for output image
    6161
    6262    psRegion imageRegion = {0, 0, 0, 0};
    63     psFits *reference = psFitsAlloc(config->refFile);
     63    psFits *reference = psFitsOpen(config->refFile, "r");
    6464    if (config->reverse) {
    65         header = psFitsReadHeader(NULL, reference);
     65        header = psFitsReadHeader(NULL, reference);
    6666    }
    6767    psImage *refImage = psFitsReadImage(NULL, reference, imageRegion, 0);
    6868    if (refImage == NULL) {
    69         psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->refFile);
    70         exit(EXIT_FAILURE);
     69        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->refFile);
     70        exit(EXIT_FAILURE);
    7171    }
    7272    psTrace("pois", 3, "Read %s: %dx%d\n", config->refFile, refImage->numCols, refImage->numRows);
    73     psFree(reference);
     73    psFitsClose(reference);
    7474    // Convert to 32-bit floating point, in necessary
    7575    if (refImage->type.type != PS_TYPE_F32) {
    76         psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->refFile);
    77         psImage *temp = psImageCopy(NULL, refImage, PS_TYPE_F32);
    78         psFree(refImage);
    79         refImage = temp;
    80     }
    81 
    82     psFits *input = psFitsAlloc(config->inFile);
     76        psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->refFile);
     77        psImage *temp = psImageCopy(NULL, refImage, PS_TYPE_F32);
     78        psFree(refImage);
     79        refImage = temp;
     80    }
     81
     82    psFits *input = psFitsOpen(config->inFile, "r");
    8383    if (! config->reverse) {
    84         header = psFitsReadHeader(NULL, input);
     84        header = psFitsReadHeader(NULL, input);
    8585    }
    8686    psImage *inImage = psFitsReadImage(NULL, input, imageRegion, 0);
    8787    if (inImage == NULL) {
    88         psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->inFile);
    89         exit(EXIT_FAILURE);
     88        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->inFile);
     89        exit(EXIT_FAILURE);
    9090    }
    9191    psTrace("pois", 3, "Read %s: %dx%d\n", config->inFile, inImage->numCols, inImage->numRows);
    92     psFree(input);
     92    psFitsClose(input);
    9393    // Convert to 32-bit floating point, in necessary
    9494    if (inImage->type.type != PS_TYPE_F32) {
    95         psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->inFile);
    96         psImage *temp = psImageCopy(NULL, inImage, PS_TYPE_F32);
    97         psFree(inImage);
    98         inImage = temp;
     95        psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->inFile);
     96        psImage *temp = psImageCopy(NULL, inImage, PS_TYPE_F32);
     97        psFree(inImage);
     98        inImage = temp;
    9999    }
    100100
    101101    if ((refImage->numCols != inImage->numCols) || (refImage->numRows != inImage->numRows)) {
    102         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    103                 "Reference and input images must have the same dimensions: %dx%d vs %dx%d\n",
    104                 refImage->numCols, refImage->numRows, inImage->numCols, inImage->numRows);
    105         exit(EXIT_FAILURE);
     102        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     103                "Reference and input images must have the same dimensions: %dx%d vs %dx%d\n",
     104                refImage->numCols, refImage->numRows, inImage->numCols, inImage->numRows);
     105        exit(EXIT_FAILURE);
    106106    }
    107107    config->xImage = refImage->numCols;
    … …  
    115115    // as the square root of the value.
    116116    if (config->background != 0.0) {
    117         (void)psBinaryOp(refImage, refImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
    118         (void)psBinaryOp(inImage, inImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
     117        (void)psBinaryOp(refImage, refImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
     118        (void)psBinaryOp(inImage, inImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
    119119    }
    120120
    … …  
    127127#ifdef TESTING
    128128    // Write the mask out
    129     psFits *maskFile = psFitsAlloc("mask.fits");
     129    psFits *maskFile = psFitsOpen("mask.fits", "w");
    130130    if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
    131         psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n");
     131        psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n");
    132132    }
    133133    psTrace("pois", 1, "Mask written to mask.fits\n");
    134     psFree(maskFile);
    135 #endif
    136 
    137     psArray *stamps = NULL;             // Array of stamps
    138     psVector *solution = NULL;          // Solution vector
    139 
    140     FILE *stampsFP = NULL;              // File pointer for stamps
     134    psFitsClose(maskFile);
     135#endif
     136
     137    psArray *stamps = NULL;             // Array of stamps
     138    psVector *solution = NULL;          // Solution vector
     139
     140    FILE *stampsFP = NULL;              // File pointer for stamps
    141141    if (config->stampFile) {
    142         psTrace("pois", 5, "Opening stamp file, %s\n", config->stampFile);
    143         stampsFP = fopen(config->stampFile, "r");
    144         if (! stampsFP) {
    145             psError(PS_ERR_IO, true, "Unable to open stamps file, %s!\n", config->stampFile);
    146             exit(EXIT_FAILURE);
    147         }
     142        psTrace("pois", 5, "Opening stamp file, %s\n", config->stampFile);
     143        stampsFP = fopen(config->stampFile, "r");
     144        if (! stampsFP) {
     145            psError(PS_ERR_IO, true, "Unable to open stamps file, %s!\n", config->stampFile);
     146            exit(EXIT_FAILURE);
     147        }
    148148    }
    149149
    150150    // Iterate for a good solution
    151     bool badStamps = true;              // Do we have bad stamps, such that we need to continue to iterate?
     151    bool badStamps = true;              // Do we have bad stamps, such that we need to continue to iterate?
    152152    for (int iterNum = 0; iterNum < config->numIter && badStamps; iterNum++) {
    153         psTrace("pois", 1, "Iteration %d...\n", iterNum);
    154 
    155         // Find stamps
    156         if (config->stampFile) {
    157             stamps = poisReadStamps(stamps, &stampsFP, refImage, mask, config);
    158         } else {
    159             stamps = poisFindStamps(stamps, refImage, mask, config);
    160         }
    161         int numStamps = 0;
    162         for (int s = 0; s < stamps->n; s++) {
    163             poisStamp *stamp = stamps->data[s]; // Stamp of interest
    164             if (stamp->status != POIS_STAMP_BAD) {
    165                 numStamps++;
    166             }
    167         }
    168         psTrace("pois.time", 1, "%d stamps found at %f sec\n", numStamps, getTime() - startTime);
    169         if (numStamps == 0) {
    170             fprintf(stderr, "No stamps found.  Check input parameters.\n");
    171             exit(EXIT_FAILURE);
    172         }
    173 
    174         // Calculate equation
    175         (void)poisCalculateEquation(stamps, refImage, inImage, kernelBasisFunctions, config);
    176         psTrace("pois.time", 1, "Equation calculated at %f sec\n", getTime() - startTime);
    177 
    178         // Solve the equation
    179         solution = poisSolveEquation(solution, stamps, config);
    180         psTrace("pois.time", 1, "Matrix equation solved at %f sec\n", getTime() - startTime);
    181 
    182 #ifdef TESTING
    183         // Print solution
    184         psTrace("pois.solution", 5, "Solution:\n");
    185         for (int i = 0; i < solution->n - 1; i++) {
    186             poisKernelBasis *kernel = kernelBasisFunctions->data[i]; // The i-th kernel basis function
    187             psTrace("pois.solution", 7, "%d: (%d,%d) x^%d y^%d --> %f\n", i, kernel->u, kernel->v,
    188                     kernel->xOrder, kernel->yOrder, solution->data.F64[i]);
    189         }
    190 #endif
    191 
    192 #ifdef TESTING
    193         // Save the kernel postage stamp
    194         char kernelName[MAXCHAR];               // File name for kernel
    195         snprintf(kernelName, MAXCHAR, "kernel%d.fits", iterNum);
    196         psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config);
    197         psFits *kernelFile = psFitsAlloc(kernelName);
    198         if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0)) {
    199             psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName);
    200         }
    201         psTrace("pois", 1, "Kernel written to %s\n", kernelName);
    202         psFree(kernelFile);
    203         psFree(kernelImage);
    204 #endif
    205 
    206         // Calculate deviations for the stamps
    207         psVector *deviations = poisCalculateDeviations(NULL, stamps, refImage, inImage, mask,
    208                                                        kernelBasisFunctions, solution, config);
    209         badStamps = poisRejectStamps(stamps, mask, deviations, config);
    210 
    211         psFree(deviations);
     153        psTrace("pois", 1, "Iteration %d...\n", iterNum);
     154
     155        // Find stamps
     156        if (config->stampFile) {
     157            stamps = poisReadStamps(stamps, &stampsFP, refImage, mask, config);
     158        } else {
     159            stamps = poisFindStamps(stamps, refImage, mask, config);
     160        }
     161        int numStamps = 0;
     162        for (int s = 0; s < stamps->n; s++) {
     163            poisStamp *stamp = stamps->data[s]; // Stamp of interest
     164            if (stamp->status != POIS_STAMP_BAD) {
     165                numStamps++;
     166            }
     167        }
     168        psTrace("pois.time", 1, "%d stamps found at %f sec\n", numStamps, getTime() - startTime);
     169        if (numStamps == 0) {
     170            fprintf(stderr, "No stamps found.  Check input parameters.\n");
     171            exit(EXIT_FAILURE);
     172        }
     173
     174        // Calculate equation
     175        (void)poisCalculateEquation(stamps, refImage, inImage, kernelBasisFunctions, config);
     176        psTrace("pois.time", 1, "Equation calculated at %f sec\n", getTime() - startTime);
     177
     178        // Solve the equation
     179        solution = poisSolveEquation(solution, stamps, config);
     180        psTrace("pois.time", 1, "Matrix equation solved at %f sec\n", getTime() - startTime);
     181
     182#ifdef TESTING
     183        // Print solution
     184        psTrace("pois.solution", 5, "Solution:\n");
     185        for (int i = 0; i < solution->n - 1; i++) {
     186            poisKernelBasis *kernel = kernelBasisFunctions->data[i]; // The i-th kernel basis function
     187            psTrace("pois.solution", 7, "%d: (%d,%d) x^%d y^%d --> %f\n", i, kernel->u, kernel->v,
     188                    kernel->xOrder, kernel->yOrder, solution->data.F64[i]);
     189        }
     190#endif
     191
     192#ifdef TESTING
     193        // Save the kernel postage stamp
     194        char kernelName[MAXCHAR];               // File name for kernel
     195        snprintf(kernelName, MAXCHAR, "kernel%d.fits", iterNum);
     196        psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config);
     197        psFits *kernelFile = psFitsOpen(kernelName, "w");
     198        if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0)) {
     199            psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName);
     200        }
     201        psTrace("pois", 1, "Kernel written to %s\n", kernelName);
     202        psFitsClose(kernelFile);
     203        psFree(kernelImage);
     204#endif
     205
     206        // Calculate deviations for the stamps
     207        psVector *deviations = poisCalculateDeviations(NULL, stamps, refImage, inImage, mask,
     208                                                       kernelBasisFunctions, solution, config);
     209        badStamps = poisRejectStamps(stamps, mask, deviations, config);
     210
     211        psFree(deviations);
    212212    } // Iterating for good solution
    213213
    214214    // If there was rejection on the last round, re-solve the equation using only the good stamps
    215215    if (badStamps) {
    216         solution = poisSolveEquation(solution, stamps, config);
     216        solution = poisSolveEquation(solution, stamps, config);
    217217    }
    218218
    219219    if (stampsFP) {
    220         fclose(stampsFP);
     220        fclose(stampsFP);
    221221    }
    222222
    223223    // Undo the mask for bad stamp area
    224224    for (int j = 0; j < mask->numRows; j++) {
    225         for (int i = 0; i < mask->numCols; i++) {
    226             if (mask->data.U8[j][i] & POIS_MASK_STAMP) {
    227                 mask->data.U8[j][i] &= ~POIS_MASK_STAMP;
    228             }
    229         }
     225        for (int i = 0; i < mask->numCols; i++) {
     226            if (mask->data.U8[j][i] & POIS_MASK_STAMP) {
     227                mask->data.U8[j][i] &= ~POIS_MASK_STAMP;
     228            }
     229        }
    230230    }
    231231
    … …  
    239239#ifdef TESTING
    240240    // Save the convolved image
    241     char convName[MAXCHAR];             // Filename for convolved image
     241    char convName[MAXCHAR];             // Filename for convolved image
    242242    snprintf(convName, MAXCHAR, "%s.conv", config->outFile);
    243     psFits *convFile = psFitsAlloc(convName);
     243    psFits *convFile = psFitsOpen(convName, "w");
    244244    if (!psFitsWriteImage(convFile, NULL, convImage, 0)) {
    245         psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName);
    246     }
    247     psFree(convFile);
     245        psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName);
     246    }
     247    psFitsClose(convFile);
    248248    psTrace("pois", 1, "Convolved image written to %s\n", convName);
    249249#endif
    … …  
    252252    psImage *subImage = psImageAlloc(config->xImage, config->yImage, PS_TYPE_F32);
    253253    if (config->reverse) {
    254         (void)psBinaryOp(subImage, convImage, "-", inImage);
     254        (void)psBinaryOp(subImage, convImage, "-", inImage);
    255255    } else {
    256         (void)psBinaryOp(subImage, inImage, "-", convImage);
     256        (void)psBinaryOp(subImage, inImage, "-", convImage);
    257257    }
    258258    psTrace("pois.time", 1, "Subtraction completed at %f sec\n", getTime() - startTime);
    … …  
    260260#ifdef TESTING
    261261    {
    262         psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    263         (void)psImageStats(stats, subImage, mask, 0);
    264         psTrace("pois", 5, "Subtracted image statistics: %f %f\n", stats->sampleMean, stats->sampleStdev);
    265         psFree(stats);
     262        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     263        (void)psImageStats(stats, subImage, mask, 0);
     264        psTrace("pois", 5, "Subtracted image statistics: %f %f\n", stats->sampleMean, stats->sampleStdev);
     265        psFree(stats);
    266266    }
    267267#endif
    … …  
    269269    // Mask out the edges
    270270    for (int y = config->yKernel; y < subImage->numRows - config->yKernel; y++) {
    271         for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
    272 
    273             for (int x = 0; x < config->xKernel; x++) {
    274                 subImage->data.F32[y][x] = 0.0;
    275             }
    276             for (int x = subImage->numCols - config->xKernel; x < subImage->numCols; x++) {
    277                 subImage->data.F32[y][x] = 0.0;
    278             }
    279         }
     271        for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
     272
     273            for (int x = 0; x < config->xKernel; x++) {
     274                subImage->data.F32[y][x] = 0.0;
     275            }
     276            for (int x = subImage->numCols - config->xKernel; x < subImage->numCols; x++) {
     277                subImage->data.F32[y][x] = 0.0;
     278            }
     279        }
    280280    }
    281281    for (int y = 0; y < config->yKernel; y++) {
    282         for (int x = 0; x < subImage->numCols; x++) {
    283             subImage->data.F32[y][x] = 0.0;
    284         }
     282        for (int x = 0; x < subImage->numCols; x++) {
     283            subImage->data.F32[y][x] = 0.0;
     284        }
    285285    }
    286286    for (int y = subImage->numRows - config->yKernel; y < subImage->numRows; y++) {
    287         for (int x = 0; x < subImage->numCols; x++) {
    288             subImage->data.F32[y][x] = 0.0;
    289         }
     287        for (int x = 0; x < subImage->numCols; x++) {
     288            subImage->data.F32[y][x] = 0.0;
     289        }
    290290    }
    291291
    … …  
    293293    // Mask out bad pixels
    294294    for (int y = config->yKernel; y < subImage->numRows - config->yKernel; y++) {
    295         for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
    296             if (mask->data.U8[y][x]) {
    297                 subImage->data.F32[y][x] = 0.0;
    298             }
    299         }
     295        for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
     296            if (mask->data.U8[y][x]) {
     297                subImage->data.F32[y][x] = 0.0;
     298            }
     299        }
    300300    }
    301301
    302302    if (config->mask) {
    303         char maskName[80];              // Name of mask image
    304         sprintf(maskName, "%s.mask", config->outFile);
    305         psFits *maskFile = psFitsAlloc(maskName);
     303        char maskName[80];              // Name of mask image
     304        sprintf(maskName, "%s.mask", config->outFile);
     305        psFits *maskFile = psFitsOpen(maskName, "w");
    306306        if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
    307307            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
    308308        }
    309309        psTrace("pois", 1, "Mask image written to %s\n", maskName);
    310         psFree(maskFile);
     310        psFitsClose(maskFile);
    311311    }
    312312
    … …  
    314314    int numNaN = psImageClipNaN(subImage, 0.0);
    315315    if (numNaN > 0) {
    316         printf("Masked %d NAN pixels to zero.\n", numNaN);
    317     }
    318 
    319     psFits *subFile = psFitsAlloc(config->outFile);
     316        printf("Masked %d NAN pixels to zero.\n", numNaN);
     317    }
     318
     319    psFits *subFile = psFitsOpen(config->outFile, "w");
    320320    if (!psFitsWriteImage(subFile, header, subImage, 0)) {
    321         psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile);
    322     }
    323     psFree(subFile);
     321        psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile);
     322    }
     323    psFitsClose(subFile);
    324324    psTrace("pois", 1, "Subtracted image written to %s\n", config->outFile);
    325325    psTrace("pois.time", 1, "I/O completed at %f sec\n", getTime() - startTime);

  • trunk/pois/src/poisCalculateDeviations.c

    r5717 r5742  
    66#define MAXCHAR 80
    77
    8 psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL
    9                                   psArray *stamps, // Array of stamps
    10                                   psImage *refImage, // Reference image
    11                                   psImage *inImage, // Input image
    12                                   psImage *mask, // Mask image
    13                                   psArray *kernelParams, // Array of kernel parameters
    14                                   psVector *solution, // Solution vector
    15                                   poisConfig *config // Configuration
     8psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL
     9                                  psArray *stamps, // Array of stamps
     10                                  psImage *refImage, // Reference image
     11                                  psImage *inImage, // Input image
     12                                  psImage *mask, // Mask image
     13                                  psArray *kernelParams, // Array of kernel parameters
     14                                  psVector *solution, // Solution vector
     15                                  poisConfig *config // Configuration
    1616    )
    1717{
    … …  
    3232
    3333    if (!deviations) {
    34         deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);
     34        deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);
    3535    }
    3636
    37     int footprint = config->footprint;  // Size of stamp footprint
     37    int footprint = config->footprint;  // Size of stamp footprint
    3838    int xSize = footprint + config->xKernel; // (Half) size of subimage in x
    3939    int ySize = footprint + config->yKernel; // (Half) size of subimage in y
    40     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics
     40    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics
    4141
    4242    psImage *subStamp = psImageAlloc(2 * xSize, 2 * ySize, PS_TYPE_F32); // Subtraction of stamp
    4343    for (int s = 0; s < stamps->n; s++) {
    44         poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest
    45         int x = stamp->x;               // Stamp x coord
    46         int y = stamp->y;               // Stamp y coord
    47         if (stamp->status == POIS_STAMP_USED) {
    48             psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize};
    49             psImage *refStamp = psImageSubset(refImage, stampRegion);
    50             psImage *inStamp = psImageSubset(inImage, stampRegion);
    51             psImage *maskStamp = psImageSubset(mask, stampRegion);
    52             psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config);
    53             // Calculate chi^2
    54             (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp);
    55             (void)psBinaryOp(subStamp, subStamp, "*", subStamp);
    56             (void)psBinaryOp(subStamp, subStamp, "/", inStamp);
    57             psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,
    58                                    config->yKernel + 2 * footprint };
    59             psImage *subStampTrim = psImageSubset(subStamp, stampTrim);
    60             psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);
    61             // Copy image to workaround bug 305
    62             psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32);
    63             psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8);
    64            
    65             (void)psImageStats(stats, tempImage, tempMask, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);
    66            
    67             psFree(tempImage);
    68             psFree(tempMask);
    69            
    70             deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0);
    71             psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y,
    72                     deviations->data.F32[s]);
    73            
     44        poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest
     45        int x = stamp->x;               // Stamp x coord
     46        int y = stamp->y;               // Stamp y coord
     47        if (stamp->status == POIS_STAMP_USED) {
     48            psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize};
     49            psImage *refStamp = psImageSubset(refImage, stampRegion);
     50            psImage *inStamp = psImageSubset(inImage, stampRegion);
     51            psImage *maskStamp = psImageSubset(mask, stampRegion);
     52            psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config);
     53            // Calculate chi^2
     54            (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp);
     55            (void)psBinaryOp(subStamp, subStamp, "*", subStamp);
     56            (void)psBinaryOp(subStamp, subStamp, "/", inStamp);
     57            psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,
     58                                   config->yKernel + 2 * footprint };
     59            psImage *subStampTrim = psImageSubset(subStamp, stampTrim);
     60            psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);
     61            // Copy image to workaround bug 305
     62            psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32);
     63            psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8);
     64
     65            (void)psImageStats(stats, tempImage, tempMask, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);
     66
     67            psFree(tempImage);
     68            psFree(tempMask);
     69
     70            deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0);
     71            psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y,
     72                    deviations->data.F32[s]);
     73
    7474#ifdef TESTING
    75             char stampName[MAXCHAR];            // File name for stamp
    76             snprintf(stampName, MAXCHAR, "stamp%d.fits", s);
    77             psFits *stampFile = psFitsAlloc(stampName);
    78             if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) {
    79                 psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);
    80             }
    81             psFree(stampFile);
     75            char stampName[MAXCHAR];            // File name for stamp
     76            snprintf(stampName, MAXCHAR, "stamp%d.fits", s);
     77            psFits *stampFile = psFitsOpen(stampName, "w");
     78            if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) {
     79                psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);
     80            }
     81            psFitsClose(stampFile);
    8282#endif
    8383
    8484#if 0
    85             psFree(convRefStamp);
    86             psFree(maskStampTrim);
    87             psFree(subStampTrim);
    88             psFree(maskStamp);
    89             psFree(refStamp);
    90             psFree(inStamp);
     85            psFree(convRefStamp);
     86            psFree(maskStampTrim);
     87            psFree(subStampTrim);
     88            psFree(maskStamp);
     89            psFree(refStamp);
     90            psFree(inStamp);
    9191#endif
    92         }
     92        }
    9393    }
    9494
Note: See TracChangeset for help on using the changeset viewer.