IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5717 for trunk/pois/src/pois.c


Ignore:
Timestamp:
Dec 6, 2005, 6:43:52 PM (21 years ago)
Author:
Paul Price
Message:

Importing current working PAP version (many bug fixes since this branch) straight on top of this version, and tidying up a bit

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/pois/src/pois.c

    r3813 r5717  
    77
    88#include <stdio.h>
     9#include <string.h>
    910#include <sys/time.h>
    1011#include "pslib.h"
    … …  
    2526int main(int argc, char **argv)
    2627{
    27     poisErrorRegister();
    28 
    2928    double startTime = getTime();
    3029
    3130    // Set trace levels
    3231    psTraceSetLevel(".", 0);
    33     psTraceSetLevel("pois", 10);
    34     psTraceSetLevel("pois.config", 10);
    35     psTraceSetLevel("pois.time", 10);
    36     psTraceSetLevel("pois.solution", 0);
    37     psTraceSetLevel("pois.kernelBasisFunctions", 10);
    38     psTraceSetLevel("pois.calculateEquation", 10);
    39     psTraceSetLevel("pois.convolveImage", 10);
    40     psTraceSetLevel("pois.findStamps", 10);
    41     psTraceSetLevel("pois.parseConfig", 10);
    42     psTraceSetLevel("pois.makeMask", 10);
    43     psTraceSetLevel("pois.extractKernel", 10);
    44     psTraceSetLevel("pois.calculateDeviations", 10);
    45 
    46     // Set logging level
    47     psLogSetLevel(9);
    4832
    4933    // Parse the command line
    5034    poisConfig *config = poisParseConfig(argc, argv); // Configuration values
    5135
     36    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);
     57    }
     58
    5259    // Read inputs
    53     psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0);
     60    psMetadata *header = NULL;          // Header for output image
     61
     62    psRegion imageRegion = {0, 0, 0, 0};
    5463    psFits *reference = psFitsAlloc(config->refFile);
    55     //psMetadata *refHeader = psFitsReadHeader(NULL, reference);
    56     psImage *refImage = psFitsReadImage(NULL, reference, *imageRegion, 0);
     64    if (config->reverse) {
     65        header = psFitsReadHeader(NULL, reference);
     66    }
     67    psImage *refImage = psFitsReadImage(NULL, reference, imageRegion, 0);
    5768    if (refImage == NULL) {
    5869        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->refFile);
    … …  
    6172    psTrace("pois", 3, "Read %s: %dx%d\n", config->refFile, refImage->numCols, refImage->numRows);
    6273    psFree(reference);
     74    // Convert to 32-bit floating point, in necessary
     75    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    }
    6381
    6482    psFits *input = psFitsAlloc(config->inFile);
    65     //psMetadata *inHeader = psFitsReadHeader(NULL, input);
    66     psImage *inImage = psFitsReadImage(NULL, input, *imageRegion, 0);
     83    if (! config->reverse) {
     84        header = psFitsReadHeader(NULL, input);
     85    }
     86    psImage *inImage = psFitsReadImage(NULL, input, imageRegion, 0);
    6787    if (inImage == NULL) {
    6888        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->inFile);
    … …  
    7191    psTrace("pois", 3, "Read %s: %dx%d\n", config->inFile, inImage->numCols, inImage->numRows);
    7292    psFree(input);
    73     psFree(imageRegion);
     93    // Convert to 32-bit floating point, in necessary
     94    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;
     99    }
    74100
    75101    if ((refImage->numCols != inImage->numCols) || (refImage->numRows != inImage->numRows)) {
    … …  
    102128    // Write the mask out
    103129    psFits *maskFile = psFitsAlloc("mask.fits");
    104     if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
     130    if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
    105131        psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n");
    106132    }
    … …  
    110136
    111137    psArray *stamps = NULL;             // Array of stamps
    112     psVector *stampMask = NULL;         // Mask for stamps
    113138    psVector *solution = NULL;          // Solution vector
    114139
     140    FILE *stampsFP = NULL;              // File pointer for stamps
     141    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        }
     148    }
     149
    115150    // Iterate for a good solution
    116     psList *badStamps = NULL;           // 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?
    117152    for (int iterNum = 0; iterNum < config->numIter && badStamps; iterNum++) {
    118153        psTrace("pois", 1, "Iteration %d...\n", iterNum);
    119154
    120         psFree(badStamps);              // left over from last iteration
    121        
    122155        // Find stamps
    123         stamps = poisFindStamps(stamps, refImage, mask, config);
     156        if (config->stampFile) {
     157            stamps = poisReadStamps(stamps, &stampsFP, refImage, mask, config);
     158        } else {
     159            stamps = poisFindStamps(stamps, refImage, mask, config);
     160        }
    124161        int numStamps = 0;
    125162        for (int s = 0; s < stamps->n; s++) {
    … …  
    159196        psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config);
    160197        psFits *kernelFile = psFitsAlloc(kernelName);
    161         if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0, NULL)) {
     198        if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0)) {
    162199            psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName);
    163200        }
    … …  
    177214    // If there was rejection on the last round, re-solve the equation using only the good stamps
    178215    if (badStamps) {
    179         psFree(badStamps);
    180216        solution = poisSolveEquation(solution, stamps, config);
     217    }
     218
     219    if (stampsFP) {
     220        fclose(stampsFP);
    181221    }
    182222
    … …  
    185225        for (int i = 0; i < mask->numCols; i++) {
    186226            if (mask->data.U8[j][i] & POIS_MASK_STAMP) {
    187                 mask->data.U8[j][i] = POIS_MASK_OK;
    188             }
    189         }
    190     }
     227                mask->data.U8[j][i] &= ~POIS_MASK_STAMP;
     228            }
     229        }
     230    }
     231
     232    // Check the result
     233    (void)poisCheckKernel(solution, kernelBasisFunctions, config);
    191234
    192235    // Convolve the input image
    … …  
    199242    snprintf(convName, MAXCHAR, "%s.conv", config->outFile);
    200243    psFits *convFile = psFitsAlloc(convName);
    201     if (!psFitsWriteImage(convFile, NULL, convImage, 0, NULL)) {
     244    if (!psFitsWriteImage(convFile, NULL, convImage, 0)) {
    202245        psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName);
    203246    }
    … …  
    208251    // Do the subtraction
    209252    psImage *subImage = psImageAlloc(config->xImage, config->yImage, PS_TYPE_F32);
    210     (void)psBinaryOp(subImage, inImage, "-", convImage);
     253    if (config->reverse) {
     254        (void)psBinaryOp(subImage, convImage, "-", inImage);
     255    } else {
     256        (void)psBinaryOp(subImage, inImage, "-", convImage);
     257    }
    211258    psTrace("pois.time", 1, "Subtraction completed at %f sec\n", getTime() - startTime);
    212259
    … …  
    231278            }
    232279        }
    233        
    234         for (int y = 0; y < config->yKernel; y++) {
    235             for (int x = 0; x < subImage->numCols; x++) {
    236                 subImage->data.F32[y][x] = 0.0;
    237             }
    238         }
    239         for (int y = subImage->numRows - config->yKernel; y < subImage->numRows; y++) {
    240             for (int x = 0; x < subImage->numCols; x++) {
    241                 subImage->data.F32[y][x] = 0.0;
    242             }
    243         }
    244     }
     280    }
     281    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        }
     285    }
     286    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        }
     290    }
     291
    245292
    246293    // Mask out bad pixels
    … …  
    253300    }
    254301
     302    if (config->mask) {
     303        char maskName[80];              // Name of mask image
     304        sprintf(maskName, "%s.mask", config->outFile);
     305        psFits *maskFile = psFitsAlloc(maskName);
     306        if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
     307            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
     308        }
     309        psTrace("pois", 1, "Mask image written to %s\n", maskName);
     310        psFree(maskFile);
     311    }
     312
     313
     314    int numNaN = psImageClipNaN(subImage, 0.0);
     315    if (numNaN > 0) {
     316        printf("Masked %d NAN pixels to zero.\n", numNaN);
     317    }
     318
    255319    psFits *subFile = psFitsAlloc(config->outFile);
    256     if (!psFitsWriteImage(subFile, NULL, subImage, 0, NULL)) {
     320    if (!psFitsWriteImage(subFile, header, subImage, 0)) {
    257321        psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile);
    258322    }
    … …  
    261325    psTrace("pois.time", 1, "I/O completed at %f sec\n", getTime() - startTime);
    262326
    263 
    264327    // Clean up
    265     (void)psFree(kernelBasisFunctions);
    266     (void)psFree(refImage);
    267     (void)psFree(inImage);
    268     (void)psFree(convImage);
    269     (void)psFree(subImage);
    270     (void)psFree(config);
     328    psFree(header);
     329    psFree(kernelBasisFunctions);
     330    psFree(refImage);
     331    psFree(inImage);
     332    psFree(convImage);
     333    psFree(subImage);
     334    psFree(config);
    271335}
     336
Note: See TracChangeset for help on using the changeset viewer.