IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 13, 2005, 5:47:35 PM (20 years ago)
Author:
Paul Price
Message:

Upgrading for psLib-0.9.0

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/archive/scripts/src/phase2/papPhase2.c

    r5651 r5786  
    1717#include "pmSubtractBias.h"
    1818#include "pmChipMosaic.h"
    19 #include "pmFPAMorph.h"
     19//#include "pmFPAMorph.h"
    2020
    2121// Phase 2 needs to:
    … …  
    3939
    4040// phase2 INPUT.fits OUTPUT.fits -bias BIAS.fits -dark DARK.fits -flat FLAT.fits -chip CHIP
    41 // 
     41//
    4242// Most are self-explanatory.  "-chip" says "only work on this particular chip".
    4343
    4444
    45 #define RECIPE "PHASE2"                 // Name of the recipe to use
     45#define RECIPE "PHASE2"                 // Name of the recipe to use
    4646
    4747static psMemId memPrintAlloc(const psMemBlock *mb)
    4848{
    4949    printf("Allocated memory block %ld (%ld):\n"
    50            "\tFile %s, line %d, size %d\n"
    51            "\tPosts: %x %x %x\n",
    52            mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
    53            *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
     50           "\tFile %s, line %d, size %d\n"
     51           "\tPosts: %x %x %x\n",
     52           mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
     53           *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
    5454    return 0;
    5555}
    … …  
    5757{
    5858    printf("Freed memory block %ld (%ld):\n"
    59            "\tFile %s, line %d, size %d\n"
    60            "\tPosts: %x %x %x\n",
    61            mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
    62            *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
     59           "\tFile %s, line %d, size %d\n"
     60           "\tPosts: %x %x %x\n",
     61           mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
     62           *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
    6363    return 0;
    6464}
    … …  
    6767    psMemBlock *mb = ((psMemBlock*)ptr) - 1;
    6868    printf("Memory block %ld (%ld):\n"
    69            "\tFile %s, line %d, size %d\n"
    70            "\tPosts: %x %x %x\n",
    71            mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
    72            *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
     69           "\tFile %s, line %d, size %d\n"
     70           "\tPosts: %x %x %x\n",
     71           mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
     72           *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
    7373}
    7474
    … …  
    102102    psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", 0, "Chip number to process (if positive)", -1);
    103103    if (! psArgumentParse(arguments, &argc, argv) || argc != 3) {
    104         printf("\nPan-STARRS Phase 2 processing\n\n");
    105         printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]);
    106         psArgumentHelp(arguments);
    107         psFree(arguments);
    108         exit(EXIT_FAILURE);
    109     }
    110     const char *inputName = argv[1];    // Name of input image
    111     const char *outputName = argv[2];   // Name of output image
     104        printf("\nPan-STARRS Phase 2 processing\n\n");
     105        printf("Usage: %s INPUT.fits OUTPUT.fits\n\n", argv[0]);
     106        psArgumentHelp(arguments);
     107        psFree(arguments);
     108        exit(EXIT_FAILURE);
     109    }
     110    const char *inputName = argv[1];    // Name of input image
     111    const char *outputName = argv[2];   // Name of output image
    112112    const char *biasName = psMetadataLookupString(NULL, arguments, "-bias"); // Name of bias image
    113113    const char *darkName = psMetadataLookupString(NULL, arguments, "-dark"); // Name of dark image
    … …  
    123123
    124124    // Open the input
    125 #ifdef PRODUCTION
    126125    psFits *inputFile = psFitsOpen(inputName, "r"); // File handle for FITS file
    127 #else
    128     psFits *inputFile = psFitsAlloc(inputName); // File handle for FITS file
    129 #endif
    130126    if (! inputFile) {
    131         psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName);
     127        psErrorStackPrint(stderr, "Can't open input image: %s\n", inputName);
    132128        exit(EXIT_FAILURE);
    133129    }
    134130    psMetadata *header = psFitsReadHeader(NULL, inputFile); // FITS header
    135131#if PRODUCTION
    136     psDB *database = pmConfigDB(site);  // Database handle
     132    psDB *database = pmConfigDB(site);  // Database handle
    137133#else
    138     psDB *database = NULL;              // Database handle
     134    psDB *database = NULL;              // Database handle
    139135#endif
    140136
    141137    // Open the output and output mask
    142138    // We do it here so that we don't process the whole lot and then find out we can't open the output file
    143 #ifdef PRODUCTION
    144139    psFits *outputFile = psFitsOpen(outputName, "w");
    145 #else
    146     psFits *outputFile = psFitsAlloc(outputName);
    147 #endif
    148140    if (! outputFile) {
    149         psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName);
    150         exit(EXIT_FAILURE);
     141        psErrorStackPrint(stderr, "Can't open output image: %s\n", outputName);
     142        exit(EXIT_FAILURE);
    151143    }
    152144
    153145    // Get camera configuration from header if not already defined
    154146    if (! camera) {
    155         camera = pmConfigCameraFromHeader(site, header);
    156         if (! camera) {
    157             psErrorStackPrint(stderr, "Can't find camera configuration!\n");
    158             exit(EXIT_FAILURE);
    159         }
     147        camera = pmConfigCameraFromHeader(site, header);
     148        if (! camera) {
     149            psErrorStackPrint(stderr, "Can't find camera configuration!\n");
     150            exit(EXIT_FAILURE);
     151        }
    160152   } else if (! pmConfigValidateCamera(camera, header)) {
    161         psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName);
    162         exit(EXIT_FAILURE);
     153        psError(PS_ERR_IO, true, "%s does not seem to be from the camera.\n", inputName);
     154        exit(EXIT_FAILURE);
    163155    }
    164156    if (! recipe && !(recipe = pmConfigRecipeFromCamera(camera, RECIPE))) {
    … …  
    175167
    176168    // Set various tasks
    177     bool doMask = false;                // Mask bad pixles
    178     bool doNonLin = false;              // Non-linearity correction
    179     bool doBias = false;                // Bias subtraction
    180     bool doDark = false;                // Dark subtraction
    181     bool doOverscan = false;            // Overscan subtraction
    182     bool doFlat = false;                // Flat-field normalisation
    183     bool doFringe = false;              // Fringe subtraction
    184     bool doSource = false;              // Source identification and photometry
    185     bool doAstrom = false;              // Astrometry
     169    bool doMask = false;                // Mask bad pixles
     170    bool doNonLin = false;              // Non-linearity correction
     171    bool doBias = false;                // Bias subtraction
     172    bool doDark = false;                // Dark subtraction
     173    bool doOverscan = false;            // Overscan subtraction
     174    bool doFlat = false;                // Flat-field normalisation
     175    bool doFringe = false;              // Fringe subtraction
     176    bool doSource = false;              // Source identification and photometry
     177    bool doAstrom = false;              // Astrometry
    186178    pmOverscanAxis overscanMode = PM_OVERSCAN_NONE; // Axis for overscan
    187179    pmFit overscanFitType = PM_FIT_NONE; // Fit type for overscan
    188     int overscanBins = 1;               // Number of pixels per bin for overscan
    189     psStats *overscanStats = NULL;      // Statistics for overscan
    190     void *overscanFit = NULL;           // Overscan fit (polynomial or spline)
     180    int overscanBins = 1;               // Number of pixels per bin for overscan
     181    psStats *overscanStats = NULL;      // Statistics for overscan
     182    void *overscanFit = NULL;           // Overscan fit (polynomial or spline)
    191183
    192184    if (psMetadataLookupBool(NULL, recipe, "MASK")) {
    193         if (strlen(maskName) > 0) {
    194             doMask = true;
    195         } else {
    196             psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking "
    197                      "will be performed.\n");
    198         }
     185        if (strlen(maskName) > 0) {
     186            doMask = true;
     187        } else {
     188            psLogMsg("phase2", PS_LOG_WARN, "Masking is desired, but no mask was supplied --- no masking "
     189                     "will be performed.\n");
     190        }
    199191    }
    200192    if (psMetadataLookupBool(NULL, recipe, "NONLIN")) {
    201         doNonLin = true;
     193        doNonLin = true;
    202194    }
    203195    if (psMetadataLookupBool(NULL, recipe, "BIAS")) {
    204         if (strlen(biasName) > 0) {
    205             doBias = true;
    206         } else {
    207             psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- "
    208                      "no bias subtraction will be performed.\n");
    209         }
     196        if (strlen(biasName) > 0) {
     197            doBias = true;
     198        } else {
     199            psLogMsg("phase2", PS_LOG_WARN, "Bias subtraction is desired, but no bias was supplied --- "
     200                     "no bias subtraction will be performed.\n");
     201        }
    210202    }
    211203    if (psMetadataLookupBool(NULL, recipe, "DARK")) {
    212         if (strlen(darkName) > 0) {
    213             doDark = true;
    214         } else {
    215             psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- "
    216                      "no dark subtraction will be performed.\n");
    217         }
     204        if (strlen(darkName) > 0) {
     205            doDark = true;
     206        } else {
     207            psLogMsg("phase2", PS_LOG_WARN, "Dark subtraction is desired, but no dark was supplied --- "
     208                     "no dark subtraction will be performed.\n");
     209        }
    218210    }
    219211    if (psMetadataLookupBool(NULL, recipe, "OVERSCAN")) {
    220         doOverscan = true;
    221         psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE");
    222         if (strcasecmp(mode, "INDIVIDUAL") == 0) {
    223             overscanMode = PM_OVERSCAN_ROWS;
    224             // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip
    225         } else if (strcasecmp(mode, "ALL") == 0) {
    226             overscanMode = PM_OVERSCAN_ALL;
    227         } else if (strcasecmp(mode, "NONE") != 0) {
    228             psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:"
    229                      " assuming NONE.\n", mode);
    230         }
    231         psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT");
    232         if (strcasecmp(fit, "POLYNOMIAL") == 0) {
    233             overscanFitType = PM_FIT_POLYNOMIAL;
    234             int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
    235             overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
    236         } else if (strcasecmp(fit, "SPLINE") == 0) {
    237             overscanFitType = PM_FIT_SPLINE;
    238             int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
    239             overscanFit = NULL;
    240 //          overscanFit = psSpline1DAlloc();
    241         } else if (strcasecmp(fit, "NONE") != 0) {
    242             psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:"
    243                      " assuming NONE.\n", fit);
    244         }
    245         overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN");
    246         if (overscanBins <= 0) {
    247             psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins);
    248             overscanBins = 1;
    249         }
    250         psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT");
    251         if (strcasecmp(stat, "MEAN")) {
    252             overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    253         } else if (strcasecmp(stat, "MEDIAN")) {
    254             overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    255         } else {
    256             psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat);
    257             overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);     
    258         }
     212        doOverscan = true;
     213        psString mode = psMetadataLookupString(NULL, recipe, "OVERSCAN.MODE");
     214        if (strcasecmp(mode, "INDIVIDUAL") == 0) {
     215            overscanMode = PM_OVERSCAN_ROWS;
     216            // By "ROWS", I mean either rows or columns --- will check the READDIR for each chip
     217        } else if (strcasecmp(mode, "ALL") == 0) {
     218            overscanMode = PM_OVERSCAN_ALL;
     219        } else if (strcasecmp(mode, "NONE") != 0) {
     220            psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.MODE (%s) is not one of NONE, INDIVIDUAL, or ALL:"
     221                     " assuming NONE.\n", mode);
     222        }
     223        psString fit = psMetadataLookupString(NULL, recipe, "OVERSCAN.FIT");
     224        if (strcasecmp(fit, "POLYNOMIAL") == 0) {
     225            overscanFitType = PM_FIT_POLYNOMIAL;
     226            int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
     227            overscanFit = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
     228        } else if (strcasecmp(fit, "SPLINE") == 0) {
     229            overscanFitType = PM_FIT_SPLINE;
     230            int order = psMetadataLookupS32(NULL, recipe, "OVERSCAN.ORDER"); // Order of polynomial fit
     231            overscanFit = NULL;
     232//          overscanFit = psSpline1DAlloc();
     233        } else if (strcasecmp(fit, "NONE") != 0) {
     234            psLogMsg("phase2", PS_LOG_WARN, "OVERSCAN.FIT (%s) is not one of NONE, POLYNOMIAL, or SPLINE:"
     235                     " assuming NONE.\n", fit);
     236        }
     237        overscanBins = psMetadataLookupS32(NULL, recipe, "OVERSCAN.BIN");
     238        if (overscanBins <= 0) {
     239            psErrorStackPrint(stderr, "OVERSCAN.BIN (%d) is non-positive --- assuming 1.\n", overscanBins);
     240            overscanBins = 1;
     241        }
     242        psString stat = psMetadataLookupString(NULL, recipe, "OVERSCAN.STAT");
     243        if (strcasecmp(stat, "MEAN")) {
     244            overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
     245        } else if (strcasecmp(stat, "MEDIAN")) {
     246            overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     247        } else {
     248            psErrorStackPrint(stderr, "OVERSCAN.STAT (%s) is not one of MEAN, MEDIAN: assuming MEAN\n", stat);
     249            overscanStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
     250        }
    259251    }
    260252
    261253    if (psMetadataLookupBool(NULL, recipe, "FLAT")) {
    262         if (strlen(flatName) > 0) {
    263             doFlat = true;
    264         } else {
    265             psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- "
    266                      "no flat-fielding will be performed.\n");
    267         }
     254        if (strlen(flatName) > 0) {
     255            doFlat = true;
     256        } else {
     257            psLogMsg("phase2", PS_LOG_WARN, "Flat-fielding is desired, but no flat was supplied --- "
     258                     "no flat-fielding will be performed.\n");
     259        }
    268260    }
    269261
    270262    // Chip selection
    271263    if (chipNum >= 0) {
    272         if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) ||
    273             ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) ||
    274             ! pmFPASelectChip(mask, chipNum)) {
    275             psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum);
    276             exit(EXIT_FAILURE);
    277         }
    278         psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum);
     264        if (! pmFPASelectChip(input, chipNum) || ! pmFPASelectChip(bias, chipNum) ||
     265            ! pmFPASelectChip(dark, chipNum) || ! pmFPASelectChip(flat, chipNum) ||
     266            ! pmFPASelectChip(mask, chipNum)) {
     267            psErrorStackPrint(stderr, "Chip number %d doesn't exist in camera.\n", chipNum);
     268            exit(EXIT_FAILURE);
     269        }
     270        psLogMsg("phase2", PS_LOG_INFO, "Operating only on chip %d\n", chipNum);
    279271    }
    280272
    … …  
    284276    psLogMsg("phase2", PS_LOG_INFO, "Opening input image: %s\n", inputName);
    285277    if (! pmFPARead(input, inputFile, header, database)) {
    286         psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName);
    287         exit(EXIT_FAILURE);
    288     }
    289 
    290     pmFPAPrint(input);
    291     exit(1);
    292 
    293 #ifdef PRODUCTION
     278        psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName);
     279        exit(EXIT_FAILURE);
     280    }
    294281    psFitsClose(inputFile);
    295 #else
    296     psFree(inputFile);
    297 #endif
    298282
    299283#if 0
    … …  
    302286//#else
    303287    {
    304         // Generate mask and weight frame
    305         p_pmHDU *hdu = NULL;
    306         if (input->hdu) {
    307             hdu = input->hdu;
    308         }
    309         psArray *chips = input->chips;
    310         for (int chipNum = 0; chipNum < chips->n; chipNum++) {
    311             pmChip *chip = chips->data[chipNum];
    312             if (chip->valid) {
    313                 if (chip->hdu) {
    314                     hdu = chip->hdu;
    315                 }
    316                 psArray *cells = chip->cells;
    317                 for (int cellNum = 0; cellNum < cells->n; cellNum++) {
    318                     pmCell *cell = cells->data[cellNum];
    319                     if (cell->valid) {
    320                         if (cell->hdu) {
    321                             hdu = cell->hdu;
    322                         }
    323                        
    324                         hdu->masks = psArrayAlloc(hdu->images->n);
    325                         hdu->weights = psArrayAlloc(hdu->images->n);
    326                         psArray *readouts = cell->readouts;
    327                         for (int readNum = 0; readNum < hdu->images->n; readNum++) {
    328                             psImage *image = hdu->images->data[readNum];
    329                             psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
    330                             psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
    331                             pmReadout *readout = readouts->data[readNum];
    332                             for (int j = 0; j < image->numRows; j++) {
    333                                 for (int i = 0; i < image->numCols; i++) {
    334                                     mask->data.U8[j][i] = j + i;
    335                                     weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]);
    336                                 }
    337                             }
    338                             hdu->masks->data[readNum] = mask;
    339                             hdu->weights->data[readNum] = weight;
    340                             readout->mask = psMemIncrRefCounter(mask);
    341                             readout->weight = psMemIncrRefCounter(weight);
    342                         }
    343                     }
    344                 }
    345             }
    346         }
     288        // Generate mask and weight frame
     289        p_pmHDU *hdu = NULL;
     290        if (input->hdu) {
     291            hdu = input->hdu;
     292        }
     293        psArray *chips = input->chips;
     294        for (int chipNum = 0; chipNum < chips->n; chipNum++) {
     295            pmChip *chip = chips->data[chipNum];
     296            if (chip->valid) {
     297                if (chip->hdu) {
     298                    hdu = chip->hdu;
     299                }
     300                psArray *cells = chip->cells;
     301                for (int cellNum = 0; cellNum < cells->n; cellNum++) {
     302                    pmCell *cell = cells->data[cellNum];
     303                    if (cell->valid) {
     304                        if (cell->hdu) {
     305                            hdu = cell->hdu;
     306                        }
     307
     308                        hdu->masks = psArrayAlloc(hdu->images->n);
     309                        hdu->weights = psArrayAlloc(hdu->images->n);
     310                        psArray *readouts = cell->readouts;
     311                        for (int readNum = 0; readNum < hdu->images->n; readNum++) {
     312                            psImage *image = hdu->images->data[readNum];
     313                            psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
     314                            psImage *weight = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     315                            pmReadout *readout = readouts->data[readNum];
     316                            for (int j = 0; j < image->numRows; j++) {
     317                                for (int i = 0; i < image->numCols; i++) {
     318                                    mask->data.U8[j][i] = j + i;
     319                                    weight->data.F32[j][i] = sqrtf(image->data.F32[j][i]);
     320                                }
     321                            }
     322                            hdu->masks->data[readNum] = mask;
     323                            hdu->weights->data[readNum] = weight;
     324                            readout->mask = psMemIncrRefCounter(mask);
     325                            readout->weight = psMemIncrRefCounter(weight);
     326                        }
     327                    }
     328                }
     329            }
     330        }
    347331    }
    348332#endif
    … …  
    354338    // Load the calibration frames, if required
    355339    if (doBias) {
     340        psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias
     341        psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias
     342        if (! pmConfigValidateCamera(camera, biasHeader)) {
     343            psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName);
     344            exit(EXIT_FAILURE);
     345        }
     346        psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName);
     347        if (! pmFPARead(bias, biasFile, biasHeader, database)) {
     348            psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName);
     349            exit(EXIT_FAILURE);
     350        }
     351        psFree(biasHeader);
     352        psFitsClose(biasFile);
     353    }
     354
     355    if (doDark) {
     356        psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark
     357        psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark
     358        if (! pmConfigValidateCamera(camera, darkHeader)) {
     359            psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName);
     360            exit(EXIT_FAILURE);
     361        }
     362        psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName);
     363        if (! pmFPARead(dark, darkFile, darkHeader, database)) {
     364            psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName);
     365            exit(EXIT_FAILURE);
     366        }
     367        psFree(darkHeader);
     368        psFitsClose(darkFile);
     369    }
     370
     371    if (doMask) {
     372        psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask
     373        psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask
     374        if (! pmConfigValidateCamera(camera, maskHeader)) {
     375            psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName);
     376            exit(EXIT_FAILURE);
     377        }
     378        psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName);
     379        if (! pmFPARead(mask, maskFile, maskHeader, database)) {
     380            psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName);
     381            exit(EXIT_FAILURE);
     382        }
     383        psFree(maskHeader);
     384        psFitsClose(maskFile);
     385    }
     386
     387    if (doFlat) {
     388        psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat
     389        if (! flatFile) {
     390            psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName);
     391            exit(EXIT_FAILURE);
     392        }
     393        psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat
     394        if (! pmConfigValidateCamera(camera, flatHeader)) {
     395            psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName);
     396            exit(EXIT_FAILURE);
     397        }
     398        psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName);
     399        if (! pmFPARead(flat, flatFile, flatHeader, database)) {
     400            psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName);
     401            exit(EXIT_FAILURE);
     402        }
     403        psFree(flatHeader);
     404        psFitsClose(flatFile);
     405    }
     406
     407    psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2"));
     408
     409    psArray *inputChips = input->chips; // Array of chips in input image
     410    psArray *biasChips = bias->chips;   // Array of chips in bias image
     411    psArray *darkChips = dark->chips;   // Array of chips in dark image
     412    psArray *maskChips = mask->chips;   // Array of chips in mask image
     413    psArray *flatChips = flat->chips;   // Array of chips in flat image
     414    for (int i = 0; i < inputChips->n; i++) {
     415        pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image
     416        pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image
     417        pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image
     418        pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image
     419        pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image
     420
     421        if (! inputChip->valid) {
     422            continue;
     423        }
     424
     425        psArray *inputCells = inputChip->cells; // Array of cells in input image
     426        psArray *biasCells = biasChip->cells; // Array of cells in bias image
     427        psArray *darkCells = darkChip->cells; // Array of cells in dark image
     428        psArray *maskCells = maskChip->cells; // Array of cells in mask image
     429        psArray *flatCells = flatChip->cells; // Array of cells in flat image
     430
     431        for (int j = 0; j < inputCells->n; j++) {
     432            pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image
     433            pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image
     434            pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag
     435            pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image
     436            pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image
     437
     438            if (! inputCell->valid) {
     439                continue;
     440            }
     441
     442            psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image
     443            if (doBias && biasCell->readouts->n > 1) {
     444                psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be"
     445                         " used.");
     446            }
     447            if (doDark && darkCell->readouts->n > 1) {
     448                psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be"
     449                         " used.");
     450            }
     451            if (doMask && maskCell->readouts->n > 1) {
     452                psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be"
     453                         " used.");
     454            }
     455            if (doFlat && flatCell->readouts->n > 1) {
     456                psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be"
     457                         " used.");
     458            }
     459
     460            pmReadout *biasReadout = NULL; // Bias readout
     461            pmReadout *maskReadout = NULL; // Mask readout
     462            pmReadout *flatReadout = NULL; // Flat readout
     463            psImage *biasImage = NULL;  // Bias pixels
     464            psImage *darkImage = NULL;  // Dark pixels
     465            float darkTime = 1.0;       // Dark time for dark image
     466            psImage *maskImage = NULL;  // Mask pixels
     467
     468            if (doBias) {
     469                biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image
     470                biasImage = biasReadout->image;
     471            }
     472            if (doDark) {
     473                pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image
     474                darkImage = darkReadout->image;
     475                darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME");
     476                if (darkTime <= 0.0) {
     477                    psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime);
     478                    exit(EXIT_FAILURE);
     479                }
     480            }
     481            if (doMask) {
     482                maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image
     483                maskImage = maskReadout->image;
     484            }
     485            if (doFlat) {
     486                flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image
     487            }
     488
     489            for (int k = 0; k < inputReadouts->n; k++) {
     490                pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image
     491                psImage *inputImage = inputReadout->image; // The actual image
     492
     493                // Mask bad pixels
     494                if (doMask) {
     495                    float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION");
     496
     497                    // Need to change this later to grow the mask by the size of the convolution kernel
     498                    (void)pmMaskBadPixels(inputReadout, maskImage,
     499                                          PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation,
     500                                          PM_MASK_TRAP, 0);
     501                }
     502
     503                // Non-linearity correction
     504                if (doNonLin) {
     505                    psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA");
     506                    if (! dataItem) {
     507                        psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to "
     508                                 "find NONLIN.DATA in recipe --- ignored.\n");
     509                    } else {
     510                        switch (dataItem->type) {
     511                          case PS_DATA_VECTOR:
     512                            {
     513                                // These are the polynomial coefficients
     514                                psVector *coeff = dataItem->data.V; // The coefficient vector
     515                                if (coeff->type.type != PS_TYPE_F64) {
     516                                    psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version
     517                                    coeff = temp;
     518                                }
     519                                psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1,
     520                                                                                 PS_POLYNOMIAL_ORD);
     521                                psFree(correction->coeff);
     522                                correction->coeff = psMemIncrRefCounter(coeff->data.F64);
     523                                (void)pmNonLinearityPolynomial(inputReadout, correction);
     524                                psFree(coeff);
     525                                psFree(correction);
     526                            }
     527                            break;
     528                          case PS_DATA_STRING:
     529                            {
     530                                // This is a filename
     531                                psString name = dataItem->data.V;       // Filename
     532                                psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0);
     533                                if (psLookupTableRead(table) <= 0) {
     534                                    psErrorStackPrint(stderr, "Unable to read non-linearity correction file "
     535                                                      "%s --- ignored\n", name);
     536                                } else {
    356537#ifdef PRODUCTION
    357         psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias
     538                                    (void)pmNonLinearityLookup(inputReadout, table);
    358539#else
    359         psFits *biasFile = psFitsAlloc(biasName);
    360 #endif
    361         psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias
    362         if (! pmConfigValidateCamera(camera, biasHeader)) {
    363             psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName);
    364             exit(EXIT_FAILURE);
    365         }
    366         psLogMsg("phase2", PS_LOG_INFO, "Opening bias image: %s\n", biasName);
    367         if (! pmFPARead(bias, biasFile, biasHeader, database)) {
    368             psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName);
    369             exit(EXIT_FAILURE);
    370         }
    371         psFree(biasHeader);
     540                                    psVector *influx = table->values->data[0];
     541                                    psVector *outflux = table->values->data[1];
     542                                    (void)pmNonLinearityLookup(inputReadout, table->values->data[0],
     543                                                               table->values->data[1]);
     544#endif
     545                                }
     546                                psFree(table);
     547                            }
     548                            break;
     549                          case PS_DATA_METADATA:
     550                            {
     551                                // This is a menu
     552                                psMetadata *options = dataItem->data.V; // Options with concept values as keys
     553                                bool mdok = false; // Success of MD lookup
     554                                psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE");
     555                                if (! mdok || ! concept) {
     556                                    psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
     557                                             "unable to find NONLIN.SOURCE in recipe --- ignored.\n");
     558                                } else {
     559                                    psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept);
     560                                    if (! conceptValueItem) {
     561                                        psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s "
     562                                                 "for non-linearity correction --- ignored.\n", concept);
     563                                    } else if (conceptValueItem->type != PS_DATA_STRING) {
     564                                        psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as"
     565                                                 " expected for non-linearity correction --- ignored.\n",
     566                                                 concept);
     567                                    } else {
     568                                        psString conceptValue = conceptValueItem->data.V;
     569                                        // Get the value of the concept
     570                                        psMetadataItem *optionItem = psMetadataLookup(options, conceptValue);
     571                                        if (!optionItem) {
     572                                            psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA"
     573                                                     " --- ignored.\n", conceptValue);
     574                                        } else if (optionItem->type == PS_DATA_VECTOR) {
     575                                            // These are the polynomial coefficients
     576                                            psVector *coeff = optionItem->data.V; // The coefficient vector
     577                                            if (coeff->type.type != PS_TYPE_F64) {
     578                                                psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64);
     579                                                coeff = temp;
     580                                            }
     581                                            // Polynomial correction
     582                                            psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1,
     583                                                                                             PS_POLYNOMIAL_ORD);
     584                                            psFree(correction->coeff);
     585                                            correction->coeff = psMemIncrRefCounter(coeff->data.F64);
     586                                            (void)pmNonLinearityPolynomial(inputReadout, correction);
     587                                            psFree(coeff);
     588                                            psFree(correction);
     589                                        } else if (optionItem->type == PS_DATA_STRING) {
     590                                            // This is a filename
     591                                            psString tableName = optionItem->data.V; // The filename
     592                                            psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0);
     593                                            int numLines = 0; // Number of lines read from table
     594                                            if ((numLines = psLookupTableRead(table)) <= 0) {
     595                                                psErrorStackPrint(stderr, "Unable to read non-linearity "
     596                                                                  "correction file %s --- ignored\n",
     597                                                                  tableName);
     598                                            } else {
    372599#ifdef PRODUCTION
    373         psFitsClose(biasFile);
     600                                                (void)pmNonLinearityLookup(inputReadout, table);
    374601#else
    375         psFree(biasFile);
    376 #endif
    377     }
    378 
    379     if (doDark) {
     602                                                printf("XXX: Non-linearity correction from lookup table not "
     603                                                       "yet implemented.\n");
     604#endif
     605                                            }
     606                                            psFree(table);
     607                                        } else {
     608                                            psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction "
     609                                                     "desired but unable to interpret NONLIN.DATA for %s"
     610                                                     " --- ignored\n", conceptValue);
     611                                        }
     612                                    }
     613                                }
     614                            }
     615                            break;
     616                          default:
     617                            psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
     618                                     "NONLIN.DATA is of invalid type --- ignored.\n");
     619                        }
     620                    }
     621                } // Non-linearity correction
     622
     623                // Bias, dark and overscan subtraction are all merged.
     624
     625                // Changed the definition of pmSubtractBias to do the dark subtraction as well, but
     626                // it's not in there yet.  Once it is, we can get rid of "pedestal".
     627
    380628#ifdef PRODUCTION
    381         psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark
     629                psImage *pedestal = NULL;       // Pedestal image (bias + scaled dark)
     630                if (doDark) {
     631                    // Dark time for input image
     632                    float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME");
     633                    if (inputTime <= 0) {
     634                        psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n",
     635                                          inputTime);
     636                        exit(EXIT_FAILURE);
     637                    }
     638
     639                    pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*",
     640                                                    psScalarAlloc(inputTime * darkTime, PS_TYPE_F32));
     641                }
     642                if (doBias) {
     643                    if (pedestal) {
     644                        if (biasImage->numRows != darkImage->numRows ||
     645                            biasImage->numCols != darkImage->numCols) {
     646                            psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: "
     647                                    "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows,
     648                                    darkImage->numCols, darkImage->numRows);
     649                            exit(EXIT_FAILURE);
     650                        }
     651                        (void)psBinaryOp(pedestal, pedestal, "+", biasImage);
     652                    } else {
     653                        pedestal = psMemIncrRefCounter(biasImage);
     654                    }
     655                }
     656#endif
     657                if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) {
     658                    // Need to get the read direction
     659                    int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR");
     660                    if (readdir == 1) {
     661// psmodule-0.8.0 has confused PM_OVERSCAN_ROWS and PM_OVERSCAN_COLUMNS; bug 608.
     662#ifdef PRODUCTION
     663                        overscanMode = PM_OVERSCAN_ROWS;
    382664#else
    383         psFits *darkFile = psFitsAlloc(darkName);
    384 #endif
    385         psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark
    386         if (! pmConfigValidateCamera(camera, darkHeader)) {
    387             psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName);
    388             exit(EXIT_FAILURE);
    389         }
    390         psLogMsg("phase2", PS_LOG_INFO, "Opening dark image: %s\n", darkName);
    391         if (! pmFPARead(dark, darkFile, darkHeader, database)) {
    392             psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName);
    393             exit(EXIT_FAILURE);
    394         }
    395         psFree(darkHeader);
     665                        overscanMode = PM_OVERSCAN_COLUMNS;
     666#endif
     667                    } else if (readdir == 2) {
    396668#ifdef PRODUCTION
    397         psFitsClose(darkFile);
     669                        overscanMode = PM_OVERSCAN_COLUMNS;
    398670#else
    399         psFree(darkFile);
    400 #endif
    401     }
    402 
    403     if (doMask) {
     671                        overscanMode = PM_OVERSCAN_ROWS;
     672#endif
     673                    } else {
     674                        psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n",
     675                                          readdir);
     676                        overscanMode = PM_OVERSCAN_ROWS;
     677                    }
     678                }
     679
     680                if (doBias || doOverscan || doDark) {
     681                    psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions
    404682#ifdef PRODUCTION
    405         psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask
     683                    (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode,
     684                                         overscanStats, overscanBins, overscanFitType, pedestal);
    406685#else
    407         psFits *maskFile = psFitsAlloc(maskName);
    408 #endif
    409         psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask
    410         if (! pmConfigValidateCamera(camera, maskHeader)) {
    411             psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName);
    412             exit(EXIT_FAILURE);
    413         }
    414         psLogMsg("phase2", PS_LOG_INFO, "Opening mask image: %s\n", maskName);
    415         if (! pmFPARead(mask, maskFile, maskHeader, database)) {
    416             psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName);
    417             exit(EXIT_FAILURE);
    418         }
    419         psFree(maskHeader);
    420 
    421 #ifdef PRODUCTION
    422         psFitsClose(maskFile);
    423 #else
    424         psFree(maskFile);
    425 #endif
    426     }
    427    
    428     if (doFlat) {
    429 #ifdef PRODUCTION
    430         psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat
    431 #else
    432         psFits *flatFile = psFitsAlloc(flatName);
    433 #endif
    434         if (! flatFile) {
    435             psErrorStackPrint(stderr, "Can't open flat image: %s\n", flatName);
    436             exit(EXIT_FAILURE);
    437         }
    438         psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat
    439         if (! pmConfigValidateCamera(camera, flatHeader)) {
    440             psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName);
    441             exit(EXIT_FAILURE);
    442         }
    443         psLogMsg("phase2", PS_LOG_INFO, "Opening flat image: %s\n", flatName);
    444         if (! pmFPARead(flat, flatFile, flatHeader, database)) {
    445             psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName);
    446             exit(EXIT_FAILURE);
    447         }
    448         psFree(flatHeader);
    449 #ifdef PRODUCTION
    450         psFitsClose(flatFile);
    451 #else
    452         psFree(flatFile);
    453 #endif
    454     }
    455    
    456     psLogMsg("phase2", PS_LOG_INFO, "Input completed after %f sec.\n", psTimerMark("phase2"));
    457 
    458     psArray *inputChips = input->chips; // Array of chips in input image
    459     psArray *biasChips = bias->chips;   // Array of chips in bias image
    460     psArray *darkChips = dark->chips;   // Array of chips in dark image
    461     psArray *maskChips = mask->chips;   // Array of chips in mask image
    462     psArray *flatChips = flat->chips;   // Array of chips in flat image
    463     for (int i = 0; i < inputChips->n; i++) {
    464         pmChip *inputChip = inputChips->data[i]; // Chip of interest in input image
    465         pmChip *biasChip = biasChips->data[i]; // Chip of interest in bias image
    466         pmChip *darkChip = darkChips->data[i]; // Chip of interest in dark image
    467         pmChip *maskChip = maskChips->data[i]; // Chip of interest in mask image
    468         pmChip *flatChip = flatChips->data[i]; // Chip of interest in flat image
    469 
    470         if (! inputChip->valid) {
    471             continue;
    472         }
    473 
    474         psArray *inputCells = inputChip->cells; // Array of cells in input image
    475         psArray *biasCells = biasChip->cells; // Array of cells in bias image
    476         psArray *darkCells = darkChip->cells; // Array of cells in dark image
    477         psArray *maskCells = maskChip->cells; // Array of cells in mask image
    478         psArray *flatCells = flatChip->cells; // Array of cells in flat image
    479 
    480         for (int j = 0; j < inputCells->n; j++) {
    481             pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image
    482             pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image
    483             pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag
    484             pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image
    485             pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image
    486 
    487             if (! inputCell->valid) {
    488                 continue;
    489             }
    490 
    491             psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image
    492             if (doBias && biasCell->readouts->n > 1) {
    493                 psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be"
    494                          " used.");
    495             }
    496             if (doDark && darkCell->readouts->n > 1) {
    497                 psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be"
    498                          " used.");
    499             }
    500             if (doMask && maskCell->readouts->n > 1) {
    501                 psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be"
    502                          " used.");
    503             }
    504             if (doFlat && flatCell->readouts->n > 1) {
    505                 psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be"
    506                          " used.");
    507             }
    508 
    509             pmReadout *biasReadout = NULL; // Bias readout
    510             pmReadout *maskReadout = NULL; // Mask readout
    511             pmReadout *flatReadout = NULL; // Flat readout
    512             psImage *biasImage = NULL;  // Bias pixels
    513             psImage *darkImage = NULL;  // Dark pixels
    514             float darkTime = 1.0;       // Dark time for dark image
    515             psImage *maskImage = NULL;  // Mask pixels
    516 
    517             if (doBias) {
    518                 biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image
    519                 biasImage = biasReadout->image;
    520             }
    521             if (doDark) {
    522                 pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image
    523                 darkImage = darkReadout->image;
    524                 darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME");
    525                 if (darkTime <= 0.0) {
    526                     psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime);
    527                     exit(EXIT_FAILURE);
    528                 }
    529             }
    530             if (doMask) {
    531                 maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image
    532                 maskImage = maskReadout->image;
    533             }
    534             if (doFlat) {
    535                 flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image
    536             }
    537 
    538             for (int k = 0; k < inputReadouts->n; k++) {
    539                 pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image
    540                 psImage *inputImage = inputReadout->image; // The actual image
    541 
    542                 // Mask bad pixels
    543                 if (doMask) {
    544                     float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION");
    545 
    546                     // Need to change this later to grow the mask by the size of the convolution kernel
    547                     (void)pmMaskBadPixels(inputReadout, maskImage,
    548                                           PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation,
    549                                           PM_MASK_TRAP, 0);
    550                 }
    551 
    552                 // Non-linearity correction
    553                 if (doNonLin) {
    554                     psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA");
    555                     if (! dataItem) {
    556                         psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but unable to "
    557                                  "find NONLIN.DATA in recipe --- ignored.\n");
    558                     } else {
    559                         switch (dataItem->type) {
    560                           case PS_DATA_VECTOR:
    561                             {
    562                                 // These are the polynomial coefficients
    563                                 psVector *coeff = dataItem->data.V; // The coefficient vector
    564                                 if (coeff->type.type != PS_TYPE_F64) {
    565                                     psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version
    566                                     coeff = temp;
    567                                 }
    568                                 psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1,
    569                                                                                  PS_POLYNOMIAL_ORD);
    570                                 psFree(correction->coeff);
    571                                 correction->coeff = psMemIncrRefCounter(coeff->data.F64);
    572                                 (void)pmNonLinearityPolynomial(inputReadout, correction);
    573                                 psFree(coeff);
    574                                 psFree(correction);
    575                             }
    576                             break;
    577                           case PS_DATA_STRING:
    578                             {
    579                                 // This is a filename
    580                                 psString name = dataItem->data.V;       // Filename
    581                                 psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0);
    582                                 if (psLookupTableRead(table) <= 0) {
    583                                     psErrorStackPrint(stderr, "Unable to read non-linearity correction file "
    584                                                       "%s --- ignored\n", name);
    585                                 } else {
    586 #ifdef PRODUCTION
    587                                     (void)pmNonLinearityLookup(inputReadout, table);
    588 #else
    589                                     psVector *influx = table->values->data[0];
    590                                     psVector *outflux = table->values->data[1];
    591                                     (void)pmNonLinearityLookup(inputReadout, table->values->data[0],
    592                                                                table->values->data[1]);
    593 #endif
    594                                 }
    595                                 psFree(table);
    596                             }
    597                             break;
    598                           case PS_DATA_METADATA:
    599                             {
    600                                 // This is a menu
    601                                 psMetadata *options = dataItem->data.V; // Options with concept values as keys
    602                                 bool mdok = false; // Success of MD lookup
    603                                 psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE");
    604                                 if (! mdok || ! concept) {
    605                                     psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
    606                                              "unable to find NONLIN.SOURCE in recipe --- ignored.\n");
    607                                 } else {
    608                                     psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept);
    609                                     if (! conceptValueItem) {
    610                                         psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s "
    611                                                  "for non-linearity correction --- ignored.\n", concept);
    612                                     } else if (conceptValueItem->type != PS_DATA_STRING) {
    613                                         psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as"
    614                                                  " expected for non-linearity correction --- ignored.\n",
    615                                                  concept);
    616                                     } else {
    617                                         psString conceptValue = conceptValueItem->data.V;
    618                                         // Get the value of the concept
    619                                         psMetadataItem *optionItem = psMetadataLookup(options, conceptValue);
    620                                         if (!optionItem) {
    621                                             psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA"
    622                                                      " --- ignored.\n", conceptValue);
    623                                         } else if (optionItem->type == PS_DATA_VECTOR) {
    624                                             // These are the polynomial coefficients
    625                                             psVector *coeff = optionItem->data.V; // The coefficient vector
    626                                             if (coeff->type.type != PS_TYPE_F64) {
    627                                                 psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64);
    628                                                 coeff = temp;
    629                                             }
    630                                             // Polynomial correction
    631                                             psPolynomial1D *correction = psPolynomial1DAlloc(coeff->n - 1,
    632                                                                                              PS_POLYNOMIAL_ORD);
    633                                             psFree(correction->coeff);
    634                                             correction->coeff = psMemIncrRefCounter(coeff->data.F64);
    635                                             (void)pmNonLinearityPolynomial(inputReadout, correction);
    636                                             psFree(coeff);
    637                                             psFree(correction);
    638                                         } else if (optionItem->type == PS_DATA_STRING) {
    639                                             // This is a filename
    640                                             psString tableName = optionItem->data.V; // The filename
    641                                             psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0);
    642                                             int numLines = 0; // Number of lines read from table
    643                                             if ((numLines = psLookupTableRead(table)) <= 0) {
    644                                                 psErrorStackPrint(stderr, "Unable to read non-linearity "
    645                                                                   "correction file %s --- ignored\n",
    646                                                                   tableName);
    647                                             } else {
    648 #ifdef PRODUCTION
    649                                                 (void)pmNonLinearityLookup(inputReadout, table);
    650 #else
    651                                                 printf("XXX: Non-linearity correction from lookup table not "
    652                                                        "yet implemented.\n");
    653 #endif
    654                                             }
    655                                             psFree(table);
    656                                         } else {
    657                                             psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction "
    658                                                      "desired but unable to interpret NONLIN.DATA for %s"
    659                                                      " --- ignored\n", conceptValue);
    660                                         }
    661                                     }
    662                                 }
    663                             }
    664                             break;
    665                           default:
    666                             psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
    667                                      "NONLIN.DATA is of invalid type --- ignored.\n");
    668                         }
    669                     }
    670                 } // Non-linearity correction
    671 
    672                 // Bias, dark and overscan subtraction are all merged.
    673 
    674                 // Changed the definition of pmSubtractBias to do the dark subtraction as well, but
    675                 // it's not in there yet.  Once it is, we can get rid of "pedestal".
    676 
    677 #ifdef PRODUCTION
    678                 psImage *pedestal = NULL;       // Pedestal image (bias + scaled dark)
    679                 if (doDark) {
    680                     // Dark time for input image
    681                     float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME");
    682                     if (inputTime <= 0) {
    683                         psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n",
    684                                           inputTime);
    685                         exit(EXIT_FAILURE);
    686                     }
    687 
    688                     pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*",
    689                                                     psScalarAlloc(inputTime * darkTime, PS_TYPE_F32));
    690                 }
    691                 if (doBias) {
    692                     if (pedestal) {
    693                         if (biasImage->numRows != darkImage->numRows ||
    694                             biasImage->numCols != darkImage->numCols) {
    695                             psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: "
    696                                     "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows,
    697                                     darkImage->numCols, darkImage->numRows);
    698                             exit(EXIT_FAILURE);
    699                         }
    700                         (void)psBinaryOp(pedestal, pedestal, "+", biasImage);
    701                     } else {
    702                         pedestal = psMemIncrRefCounter(biasImage);
    703                     }
    704                 }
    705 #endif
    706                 if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) {
    707                     // Need to get the read direction
    708                     int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR");
    709                     if (readdir == 1) {
    710                         overscanMode = PM_OVERSCAN_ROWS;
    711                     } else if (readdir == 2) {
    712                         overscanMode = PM_OVERSCAN_COLUMNS;
    713                     } else {
    714                         psErrorStackPrint(stderr, "CELL.READDIR (%d) is not 1 or 2 --- assuming 1.\n",
    715                                           readdir);
    716                         overscanMode = PM_OVERSCAN_ROWS;
    717                     }
    718                 }
    719 
    720                 if (doBias || doOverscan || doDark) {
    721                     psList *inputOverscans = pmReadoutGetBias(inputReadout); // List of overscan bias regions
    722 #ifdef PRODUCTION
    723                     (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode,
    724                                          overscanStats, overscanBins, overscanFitType, pedestal);
    725 #else
    726                     (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode,
    727                                          overscanStats, overscanBins, overscanFitType, biasReadout);
    728 #endif
    729                     psFree(inputOverscans);
    730 
    731                     // Output overscan fit results, if required
    732                     if (doOverscan) {
    733                         if (! overscanFit) {
    734                             psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n");
    735                         } else {
    736                             switch (overscanFitType) {
    737                               case PM_FIT_POLYNOMIAL:
    738                                 {
    739                                     psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial
    740                                     psString coeffs = NULL;     // String containing the coefficients
    741                                     for (int i = 0; i < poly->nX; i++) {
    742                                         psStringAppend(&coeffs, "%e ", poly->coeff[i]);
    743                                     }
    744                                     psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n",
    745                                              coeffs);
    746                                     psFree(coeffs);
    747                                 }
    748                                 break;
    749                               case PM_FIT_SPLINE:
    750                                 {
    751                                     psSpline1D *spline = (psSpline1D*)overscanFit; // The spline
    752                                     psString coeffs = NULL;     // String containing the coefficients
    753                                     for (int i = 0; i < spline->n; i++) {
    754                                         psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
    755                                         psStringAppend(&coeffs, "%d: ", i);
    756                                         for (int j = 0; j < poly->nX; j++) {
    757                                             psStringAppend(&coeffs, "%e ", poly->coeff[i]);
    758                                         }
    759                                         psStringAppend(&coeffs, "\n");
    760                                     }
    761                                     psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n",
    762                                              coeffs);
    763                                     psFree(coeffs);
    764                                 }
    765                                 break;
    766                               case PM_FIT_NONE:
    767                                 break;
    768                               default:
    769                                 psAbort(__func__, "Should never get here!!!\n");
    770                             }
    771                         }
    772                     }
    773                 }
    774 
    775                 // XXX: Temporary: until the other functions are altered to do this themselves.
    776                 // Trim, so that flat, fringe etc computations are faster.
     686                    (void)pmSubtractBias(inputReadout, overscanFit, inputOverscans, overscanMode,
     687                                         overscanStats, overscanBins, overscanFitType, biasReadout);
     688#endif
     689                    psFree(inputOverscans);
     690
     691                    // Output overscan fit results, if required
     692                    if (doOverscan) {
     693                        if (! overscanFit) {
     694                            psLogMsg("phase2", PS_LOG_WARN, "No fit generated!\n");
     695                        } else {
     696                            switch (overscanFitType) {
     697                              case PM_FIT_POLYNOMIAL:
     698                                {
     699                                    psPolynomial1D *poly = (psPolynomial1D*)overscanFit; // The polynomial
     700                                    psString coeffs = NULL;     // String containing the coefficients
     701                                    for (int i = 0; i < poly->nX; i++) {
     702                                        psStringAppend(&coeffs, "%e ", poly->coeff[i]);
     703                                    }
     704                                    psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n",
     705                                             coeffs);
     706                                    psFree(coeffs);
     707                                }
     708                                break;
     709                              case PM_FIT_SPLINE:
     710                                {
     711                                    psSpline1D *spline = (psSpline1D*)overscanFit; // The spline
     712                                    psString coeffs = NULL;     // String containing the coefficients
     713                                    for (int i = 0; i < spline->n; i++) {
     714                                        psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
     715                                        psStringAppend(&coeffs, "%d: ", i);
     716                                        for (int j = 0; j < poly->nX; j++) {
     717                                            psStringAppend(&coeffs, "%e ", poly->coeff[i]);
     718                                        }
     719                                        psStringAppend(&coeffs, "\n");
     720                                    }
     721                                    psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n",
     722                                             coeffs);
     723                                    psFree(coeffs);
     724                                }
     725                                break;
     726                              case PM_FIT_NONE:
     727                                break;
     728                              default:
     729                                psAbort(__func__, "Should never get here!!!\n");
     730                            }
     731                        }
     732                    }
     733                }
     734
     735                // XXX: Temporary: until the other functions are altered to do this themselves.
     736                // Trim, so that flat, fringe etc computations are faster.
    777737#if 0
    778738#ifndef PRODUCTION
    779                 psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC");
    780                 psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion);
    781                 psFree(inputReadout->image);
    782                 inputReadout->image = trimmed;
    783 #endif
    784 #endif
    785                 // Flat-field correction
    786                 if (doFlat) {
    787                     psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n");
     739                psRegion *trimRegion = psMetadataLookupPtr(NULL, inputCell->concepts, "CELL.TRIMSEC");
     740                psImage *trimmed = psImageSubset(inputReadout->image, *trimRegion);
     741                psFree(inputReadout->image);
     742                inputReadout->image = trimmed;
     743#endif
     744#endif
     745                // Flat-field correction
     746                if (doFlat) {
     747                    psLogMsg("phase2", PS_LOG_INFO, "Performing flat field.\n");
    788748#ifndef PRODUCTION
    789                     psImage *dummyImage = psImageAlloc(inputReadout->image->numCols,
    790                                                        inputReadout->image->numRows,
    791                                                        PS_TYPE_U8);
    792                     pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1);
    793                     (void)pmFlatField(inputReadout, dummyMask, flatReadout);
    794                     psFree(dummyMask);
    795                     psFree(dummyImage);
    796 #endif
    797                 }
    798 
    799                 // Fringe subtraction
    800                 if (doFringe) {
     749                    psImage *dummyImage = psImageAlloc(inputReadout->image->numCols,
     750                                                       inputReadout->image->numRows,
     751                                                       PS_TYPE_U8);
     752                    pmReadout *dummyMask = pmReadoutAlloc(NULL, dummyImage, NULL, 0, 0, 1, 1);
     753                    (void)pmFlatField(inputReadout, dummyMask, flatReadout);
     754                    psFree(dummyMask);
     755                    psFree(dummyImage);
     756#endif
     757                }
     758
     759                // Fringe subtraction
     760                if (doFringe) {
    801761#if 0
    802                     // Placeholder
    803                     pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask,
    804                                              psU8 maskVal, int binFactor, psStats *stats, float clipSD);
    805 #endif
    806                 }
    807 
    808                 // Source identification and photometry
    809                 if (doSource) {
    810                 }
    811 
    812                 // Astrometry
    813                 if (doAstrom) {
    814                 }
    815 
    816 
    817             } // Iterating over readouts
    818         } // Iterating over cells
    819 
    820         // Testing pmChipMosaic
    821         {
    822             psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip
    823             psFits *mosaicFile = psFitsAlloc("mosaic.fits"); // FITS file to which to write
    824             psFitsWriteImage(mosaicFile, NULL, mosaic, 0);
    825             psFree(mosaicFile);
    826             psFree(mosaic);             
    827         }
    828  
     762                    // Placeholder
     763                    pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask,
     764                                             psU8 maskVal, int binFactor, psStats *stats, float clipSD);
     765#endif
     766                }
     767
     768                // Source identification and photometry
     769                if (doSource) {
     770                }
     771
     772                // Astrometry
     773                if (doAstrom) {
     774                }
     775
     776
     777            } // Iterating over readouts
     778        } // Iterating over cells
     779
     780#if 0
     781        {
     782            // Testing pmChipMosaic
     783            psImage *mosaic = pmChipMosaic(inputChip, 1, 1); // Mosaic of chip
     784            psFits *mosaicFile = psFitsOpen("mosaic.fits", "w"); // FITS file to which to write
     785            psFitsWriteImage(mosaicFile, NULL, mosaic, 0);
     786            psFitsClose(mosaicFile);
     787            psFree(mosaic);
     788        }
     789#endif
     790
    829791    } // Iterating over chips
    830792
    … …  
    832794
    833795
     796#if 0
     797    {
     798        // Testing pmFPAMorph
     799        int nBad = 0;
     800        psMetadata *newCamera = psMetadataConfigParse(NULL, &nBad, "megacam_splice.config", true);
     801        pmFPA *newFPA = pmFPAConstruct(newCamera);
     802        pmFPAMorph(newFPA, input, true, 0, 0);
     803        //pmFPAPrint(newFPA);
     804        psFits *newFile = psFitsOpen("morph.fits", "w");
     805        pmFPAWrite(newFile, newFPA, database);
     806        psFitsClose(newFile);
     807        psFree(newFPA);
     808    }
     809#endif
     810
    834811#if 1
    835     // Testing pmFPAMorph
    836     psMetadata *newCamera = psMetadataConfigParse(NULL, NULL, "megacam_splice.config", true);
    837     pmFPA *newFPA = pmFPAConstruct(newCamera);
    838     pmFPAMorph(newFPA, input, true, 0, 0);
    839     psFits *newFile = psFitsAlloc("morph.fits");
    840     pmFPAWrite(newFile, newFPA, database);
    841     psFree(newFile);
    842     psFree(newFPA);
    843 #endif
    844 
    845 #if 0
    846812    // Write the output
    847813    pmFPAWrite(outputFile, input, database);
    … …  
    852818    psLogMsg("phase2", PS_LOG_INFO, "Output completed after %f sec.\n", psTimerMark("phase2"));
    853819
    854 #ifdef PRODUCTION
    855820    psFitsClose(outputFile);
    856 #else
    857     psFree(outputFile);
    858 #endif
    859821
    860822    psFree(arguments);
    … …  
    874836    psTimerStop();
    875837
    876 #if 1
     838#if 0
    877839    psMemCheckCorruption(true);
    878840    psMemBlock **leaks = NULL;          // List of leaks
    … …  
    881843#if 1
    882844    for (int i = 0; i < nLeaks; i++) {
    883         psMemBlock *mb = leaks[i];
    884         printf("Memory leak %ld (%ld):\n"
    885                "\tFile %s, line %d, size %d\n"
    886                "\tPosts: %x %x %x\n",
    887                mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
    888                *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
     845        psMemBlock *mb = leaks[i];
     846        printf("Memory leak %ld (%ld):\n"
     847               "\tFile %s, line %d, size %d\n"
     848               "\tPosts: %x %x %x\n",
     849               mb->id, mb->refCounter, mb->file, mb->lineno, mb->userMemorySize, mb->startblock, mb->endblock,
     850               *(void**)((int8_t *)(mb + 1) + mb->userMemorySize));
    889851    }
    890852#endif
Note: See TracChangeset for help on using the changeset viewer.