IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 22, 2005, 4:54:52 PM (21 years ago)
Author:
Paul Price
Message:

Working under pslib-0.7.0, but don't have the inner loop behaving yet.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4820 r5104  
    44#include "pslib.h"
    55#include "psAdditionals.h"
    6 
    7 #include "pmSubtractBias.h"
    86
    97#include "pmFPA.h"
    108#include "pmConfig.h"
    119#include "pmFPAConstruct.h"
     10#include "pmFPARead.h"
     11#include "pmFPAConceptsGet.h"
     12#include "pmFPAWrite.h"
     13
     14#include "pmFlatField.h"
     15#include "pmMaskBadPixels.h"
     16#include "pmNonLinear.h"
     17#include "pmSubtractBias.h"
    1218
    1319// Phase 2 needs to:
    … …  
    3743#define RECIPE "PHASE2"                 // Name of the recipe to use
    3844
    39 static psMemoryId memPrintAlloc(const psMemBlock *mb)
     45static psMemId memPrintAlloc(const psMemBlock *mb)
    4046{
    4147    printf("Allocated memory block %lld (%lld):\n"
    … …  
    4652    return 0;
    4753}
    48 static psMemoryId memPrintFree(const psMemBlock *mb)
     54static psMemId memPrintFree(const psMemBlock *mb)
    4955{
    5056    printf("Freed memory block %lld (%lld):\n"
    … …  
    6975#if 0
    7076    // Hunting memory leaks
    71     psMemAllocateCallbackSetID(3539);
    72     psMemFreeCallbackSetID(3539);
     77    psMemAllocateCallbackSetID(23289);
     78    psMemFreeCallbackSetID(23289);
    7379    psMemAllocateCallbackSet(memPrintAlloc);
    7480    psMemFreeCallbackSet(memPrintFree);
    7581#endif
    7682
    77 //    psTraceSetLevel(".", 10);
     83    psTraceSetLevel("pmConfigCameraFromHeader", 10);
    7884
    7985    // Parse the configurations
    … …  
    8894    // Parse other command-line arguments
    8995    psMetadata *arguments = psMetadataAlloc(); // The arguments, with default values
    90     psMetadataAddStr(arguments, PS_LIST_TAIL, "-bias", "Name of the bias image", "");
    91     psMetadataAddStr(arguments, PS_LIST_TAIL, "-dark", "Name of the dark image", "");
    92     psMetadataAddStr(arguments, PS_LIST_TAIL, "-flat", "Name of the flat-field image", "");
    93     psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", "Name of the mask image", "");
    94     psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", "Chip number to process (if positive)", -1);
     96    psMetadataAddStr(arguments, PS_LIST_TAIL, "-bias", 0, "Name of the bias image", "");
     97    psMetadataAddStr(arguments, PS_LIST_TAIL, "-dark", 0, "Name of the dark image", "");
     98    psMetadataAddStr(arguments, PS_LIST_TAIL, "-flat", 0, "Name of the flat-field image", "");
     99    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", 0, "Name of the mask image", "");
     100    psMetadataAddS32(arguments, PS_LIST_TAIL, "-chip", 0, "Chip number to process (if positive)", -1);
    95101    if (! psArgumentParse(arguments, &argc, argv) || argc != 3) {
    96102        printf("\nPan-STARRS Phase 2 processing\n\n");
    … …  
    113119    printf("Mask: %s\n", maskName);
    114120    printf("Chip: %d\n", chipNum);
    115     psFree(arguments);
    116121
    117122    // Open the input
    … …  
    132137#endif
    133138
    134 #if 0
    135139    // Open the output and output mask
    136140    // We do it here so that we don't process the whole lot and then find out we can't open the output file
    … …  
    144148        exit(EXIT_FAILURE);
    145149    }
     150
     151#if 0
    146152    psString outputMaskName = psStringCopy(outputName);
    147153    (void)psStringAppend(&outputMaskName, ".mask");
    … …  
    158164#endif
    159165
    160        
    161166    // Get camera configuration from header if not already defined
    162167    if (! camera) {
    … …  
    176181
    177182    // Construct camera in preparation for reading
    178     pmFPA *input = pmFPAConstruct(camera, database);
    179     pmFPA *mask = pmFPAConstruct(camera, database);
    180     pmFPA *bias = pmFPAConstruct(camera, database);
    181     pmFPA *dark = pmFPAConstruct(camera, database);
    182     pmFPA *flat = pmFPAConstruct(camera, database);
     183    pmFPA *input = pmFPAConstruct(camera);
     184    pmFPA *mask = pmFPAConstruct(camera);
     185    pmFPA *bias = pmFPAConstruct(camera);
     186    pmFPA *dark = pmFPAConstruct(camera);
     187    pmFPA *flat = pmFPAConstruct(camera);
    183188
    184189    // Set various tasks
    … …  
    198203
    199204    if (psMetadataLookupBool(NULL, recipe, "MASK")) {
    200         if (maskName) {
     205        if (strlen(maskName) > 0) {
    201206            doMask = true;
    202207        } else {
    … …  
    209214    }
    210215    if (psMetadataLookupBool(NULL, recipe, "BIAS")) {
    211         if (biasName) {
     216        if (strlen(biasName) > 0) {
    212217            doBias = true;
    213218        } else {
    … …  
    217222    }
    218223    if (psMetadataLookupBool(NULL, recipe, "DARK")) {
    219         if (darkName) {
     224        if (strlen(darkName) > 0) {
    220225            doDark = true;
    221226        } else {
    … …  
    260265        }
    261266    }
     267
    262268    if (psMetadataLookupBool(NULL, recipe, "FLAT")) {
    263         if (flatName) {
     269        if (strlen(flatName) > 0) {
    264270            doFlat = true;
    265271        } else {
    … …  
    269275    }
    270276
    271 #if 0
    272277    // Chip selection
    273278    if (chipNum >= 0) {
    … …  
    282287
    283288    // Read in the input pixels
    284     if (! pmFPARead(input, inputFile)) {
     289    if (! pmFPARead(input, inputFile, header, database)) {
    285290        psErrorStackPrint(stderr, "Unable to populate camera from FITS file: %s\n", inputName);
    286291        exit(EXIT_FAILURE);
    … …  
    292297#endif
    293298
     299    pmFPAPrint(input);
    294300
    295301    // Load the calibration frames, if required
    296     if (biasName) {
     302    if (doBias) {
     303#ifdef PRODUCTION
    297304        psFits *biasFile = psFitsOpen(biasName, "r"); // File handle for bias
     305#else
     306        psFits *biasFile = psFitsAlloc(biasName);
     307#endif
    298308        psMetadata *biasHeader = psFitsReadHeader(NULL, biasFile); // FITS header for bias
    299309        if (! pmConfigValidateCamera(camera, biasHeader)) {
    300             psError("phase2", true, "Bias (%s) does not seem to be from the same camera.\n", biasName);
     310            psError(PS_ERR_IO, true, "Bias (%s) does not seem to be from the same camera.\n", biasName);
    301311            exit(EXIT_FAILURE);
    302312        }
    303313        psFree(biasHeader);
    304         if (! pmFPARead(bias, biasFile)) {
     314        if (! pmFPARead(bias, biasFile, NULL, database)) {
    305315            psErrorStackPrint(stderr, "Unable to populate bias camera from fits FITS: %s\n", biasName);
    306316            exit(EXIT_FAILURE);
    307317        }
     318#ifdef PRODUCTION
    308319        psFitsClose(biasFile);
    309     }
    310 
    311     if (darkName) {
     320#else
     321        psFree(biasFile);
     322#endif
     323    }
     324
     325    if (doDark) {
     326#ifdef PRODUCTION
    312327        psFits *darkFile = psFitsOpen(darkName, "r"); // File handle for dark
     328#else
     329        psFits *darkFile = psFitsAlloc(darkName);
     330#endif
    313331        psMetadata *darkHeader = psFitsReadHeader(NULL, darkFile); // FITS header for dark
    314332        if (! pmConfigValidateCamera(camera, darkHeader)) {
    315             psError("phase2", true, "Dark (%s) does not seem to be from the same camera.\n", darkName);
     333            psError(PS_ERR_IO, true, "Dark (%s) does not seem to be from the same camera.\n", darkName);
    316334            exit(EXIT_FAILURE);
    317335        }
    318336        psFree(darkHeader);
    319         if (! pmFPARead(dark, darkFile)) {
     337        if (! pmFPARead(dark, darkFile, NULL, database)) {
    320338            psErrorStackPrint(stderr, "Unable to populate dark camera from fits FITS: %s\n", darkName);
    321339            exit(EXIT_FAILURE);
    322340        }
     341#ifdef PRODUCTION
    323342        psFitsClose(darkFile);
    324     }
    325 
    326     if (maskName) {
     343#else
     344        psFree(darkFile);
     345#endif
     346    }
     347
     348    if (doMask) {
     349#ifdef PRODUCTION
    327350        psFits *maskFile = psFitsOpen(maskName, "r"); // File handle for mask
     351#else
     352        psFits *maskFile = psFitsAlloc(maskName);
     353#endif
    328354        psMetadata *maskHeader = psFitsReadHeader(NULL, maskFile); // FITS header for mask
    329355        if (! pmConfigValidateCamera(camera, maskHeader)) {
    330             psError("phase2", true, "Mask (%s) does not seem to be from the same camera.\n", maskName);
     356            psError(PS_ERR_IO, true, "Mask (%s) does not seem to be from the same camera.\n", maskName);
    331357            exit(EXIT_FAILURE);
    332358        }
    333359        psFree(maskHeader);
    334         if (! pmFPARead(mask, maskFile)) {
     360        if (! pmFPARead(mask, maskFile, NULL, database)) {
    335361            psErrorStackPrint(stderr, "Unable to populate mask camera from fits FITS: %s\n", maskName);
    336362            exit(EXIT_FAILURE);
    337363        }
     364#ifdef PRODUCTION
    338365        psFitsClose(maskFile);
     366#else
     367        psFree(maskFile);
     368#endif
    339369    }
    340370   
    341     if (flatName) {
     371    if (doFlat) {
     372#ifdef PRODUCTION
    342373        psFits *flatFile = psFitsOpen(flatName, "r"); // File handle for flat
     374#else
     375        psFits *flatFile = psFitsAlloc(flatName);
     376#endif
    343377        psMetadata *flatHeader = psFitsReadHeader(NULL, flatFile); // FITS header for flat
    344378        if (! pmConfigValidateCamera(camera, flatHeader)) {
    345             psError("phase2", true, "Flat (%s) does not seem to be from the same camera.\n", flatName);
     379            psError(PS_ERR_IO, true, "Flat (%s) does not seem to be from the same camera.\n", flatName);
    346380            exit(EXIT_FAILURE);
    347381        }
    348382        psFree(flatHeader);
    349         if (! pmFPARead(flat, flatFile)) {
     383        if (! pmFPARead(flat, flatFile, NULL, database)) {
    350384            psErrorStackPrint(stderr, "Unable to populate flat camera from fits FITS: %s\n", flatName);
    351385            exit(EXIT_FAILURE);
    352386        }
     387#ifdef PRODUCTION
    353388        psFitsClose(flatFile);
     389#else
     390        psFree(flatFile);
     391#endif
    354392    }
    355393   
    356     psArray *inputChips = inputs->chips; // Array of chips in input image
     394    psArray *inputChips = input->chips; // Array of chips in input image
    357395    psArray *biasChips = bias->chips;   // Array of chips in bias image
    358396    psArray *darkChips = dark->chips;   // Array of chips in dark image
    … …  
    377415
    378416        for (int j = 0; j < inputCells->n; j++) {
    379             pmCell *inputCell = inputCells->data[i]; // Cell of interest in input image
    380             pmCell *biasCell = biasCells->data[i]; // Cell of interest in bias image
    381             pmCell *darkCell = darkCells->data[i]; // Cell of interest in dark imag
    382             pmCell *maskCell = maskCells->data[i]; // Cell of interest in mask image
    383             pmCell *flatCell = flatCells->data[i]; // Cell of interest in flat image
     417            pmCell *inputCell = inputCells->data[j]; // Cell of interest in input image
     418            pmCell *biasCell = biasCells->data[j]; // Cell of interest in bias image
     419            pmCell *darkCell = darkCells->data[j]; // Cell of interest in dark imag
     420            pmCell *maskCell = maskCells->data[j]; // Cell of interest in mask image
     421            pmCell *flatCell = flatCells->data[j]; // Cell of interest in flat image
    384422
    385423            if (! inputCell->valid) {
    … …  
    388426
    389427            psArray *inputReadouts = inputCell->readouts; // Array of readouts in input image
    390             if (biasCell->readouts->n != 1) {
     428            if (doBias && biasCell->readouts->n > 1) {
    391429                psLogMsg("phase2", PS_LOG_WARN, "Bias contains multiple readouts: only the first will be"
    392430                         " used.");
    393431            }
    394             if (darkCell->readouts->n != 1) {
     432            if (doDark && darkCell->readouts->n > 1) {
    395433                psLogMsg("phase2", PS_LOG_WARN, "Dark contains multiple readouts: only the first will be"
    396434                         " used.");
    397435            }
    398             if (maskCell->readouts->n != 1) {
     436            if (doMask && maskCell->readouts->n > 1) {
    399437                psLogMsg("phase2", PS_LOG_WARN, "Mask contains multiple readouts: only the first will be"
    400438                         " used.");
    401439            }
    402             if (flatCell->readouts->n != 1) {
     440            if (doFlat && flatCell->readouts->n > 1) {
    403441                psLogMsg("phase2", PS_LOG_WARN, "Flat contains multiple readouts: only the first will be"
    404442                         " used.");
    405443            }
    406444
     445            pmReadout *biasReadout = NULL; // Bias readout
     446            pmReadout *maskReadout = NULL; // Mask readout
     447            pmReadout *flatReadout = NULL; // Flat readout
    407448            psImage *biasImage = NULL;  // Bias pixels
    408449            psImage *darkImage = NULL;  // Dark pixels
    409450            float darkTime = 1.0;       // Dark time for dark image
    410 
    411             if (biasName) {
    412                 pmReadout *biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image
     451            psImage *maskImage = NULL;  // Mask pixels
     452
     453            if (doBias) {
     454                biasReadout = biasCell->readouts->data[0]; // Readout of interest in bias image
    413455                biasImage = biasReadout->image;
    414456            }
    415             if (darkName) {
     457            if (doDark) {
    416458                pmReadout *darkReadout = darkCell->readouts->data[0]; // Readout of interest in dark image
    417459                darkImage = darkReadout->image;
    418                 darkTime = pmReadoutGetDarkTime(darkReadout);
     460                darkTime = psMetadataLookupF32(NULL, darkCell->concepts, "CELL.DARKTIME");
    419461                if (darkTime <= 0.0) {
    420462                    psErrorStackPrint(stderr, "DARKTIME for dark image (%f) is non-positive.\n", darkTime);
    … …  
    422464                }
    423465            }
    424             if (maskName) {
    425                 pmReadout *maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image
    426             }
    427             if (flatName) {
    428                 pmReadout *flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image
     466            if (doMask) {
     467                maskReadout = maskCell->readouts->data[0]; // Readout of interest in mask image
     468                maskImage = maskReadout->image;
     469            }
     470            if (doFlat) {
     471                flatReadout = flatCell->readouts->data[0]; // Readout of interest in mask image
    429472            }
    430473
    … …  
    432475                pmReadout *inputReadout = inputReadouts->data[k]; // Readout of interest in input image
    433476                psImage *inputImage = inputReadout->image; // The actual image
    434                 psList *inputBias = inputReadout->bias; // List of overscan bias regions
     477                psList *inputOverscans = inputReadout->overscans; // List of overscan bias regions
    435478
    436479                // Mask bad pixels
    437480                if (doMask) {
     481                    float saturation = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.SATURATION");
     482
    438483                    // Need to change this later to grow the mask by the size of the convolution kernel
    439                     (void)pmMaskBadPixels(inputReadout, maskReadout,
    440                                           PS_MASK_TRAP | PS_MASK_BADCOL | PS_MASK_SAT, saturation,
    441                                           PS_MASK_TRAP, 0);
     484                    (void)pmMaskBadPixels(inputReadout, maskImage,
     485                                          PM_MASK_TRAP | PM_MASK_BADCOL | PM_MASK_SAT, saturation,
     486                                          PM_MASK_TRAP, 0);
    442487                }
    443488
    … …  
    450495                    } else {
    451496                        switch (dataItem->type) {
    452                           case PS_META_VECTOR:
    453                             // These are the polynomial coefficients
    454                             psVector *coeff = dataItem->data.V; // The coefficient vector
    455                             if (coeff->type.type != PS_TYPE_F64) {
    456                                 psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version
     497                          case PS_META_VEC:
     498                            {
     499                                // These are the polynomial coefficients
     500                                psVector *coeff = dataItem->data.V; // The coefficient vector
     501                                if (coeff->type.type != PS_TYPE_F64) {
     502                                    psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64); // F64 version
     503                                    psFree(coeff);
     504                                    coeff = temp;
     505                                }
     506                                psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,
     507                                                                                 coeff->n + 1);
     508                                for (int i = 0; i < coeff->n; i++) {
     509                                    correction->coeff[i] = coeff->data.F64[i];
     510                                }
     511                                (void)pmNonLinearityPolynomial(inputReadout, correction);
    457512                                psFree(coeff);
    458                                 coeff = temp;
     513                                psFree(correction);
    459514                            }
    460                             psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,
    461                                                                              coeff->n + 1);
    462                             for (int i = 0; i < coeff->n; i++) {
    463                                 correction->coeff[i] = coeff->data.F64[i];
    464                             }
    465                             (void)pmNonLinearityPolynomial(inputReadout, correction);
    466                             psFree(coeffCopy);
    467                             psFree(coeff);
    468                             psFree(correction);
    469515                            break;
    470516                          case PS_META_STR:
    471                             // This is a filename
    472                             psString name = dataItem->data.V;   // Filename
    473                             psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0);
    474                             if (psLookupTableRead(table) <= 0) {
    475                                 psErrorStackPrint(stderr, "Unable to read non-linearity correction file %s "
    476                                                   "--- ignored\n", tableName);
    477                             } else {
    478                                 (void)pmNonLinearityLookup(inputReadout, table);
     517                            {
     518                                // This is a filename
     519                                psString name = dataItem->data.V;       // Filename
     520                                psLookupTable *table = psLookupTableAlloc(name, "%f %f", 0);
     521                                if (psLookupTableRead(table) <= 0) {
     522                                    psErrorStackPrint(stderr, "Unable to read non-linearity correction file "
     523                                                      "%s --- ignored\n", name);
     524                                } else {
     525#ifdef PRODUCTION
     526                                    (void)pmNonLinearityLookup(inputReadout, table);
     527#else
     528                                    printf("XXX: Non-linearity with lookup table not yet implemented.\n");
     529#endif
     530                                }
     531                                psFree(table);
    479532                            }
    480                             psFree(table);
    481533                            break;
    482534                          case PS_META_META:
    483                             // This is a menu
    484                             psMetadata *options = dataItem->data.V; // Options with concept values as keys
    485                             bool mdok = false; // Success of MD lookup
    486                             psString concept = psMetadataLookupStr(&mdok, recipe, "NONLIN.SOURCE");
    487                             if (! mdok || ! concept) {
    488                                 psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
    489                                          "unable to find NONLIN.SOURCE in recipe --- ignored.\n");
    490                             } else {
    491                                 psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept);
    492                                 if (! conceptValueItem) {
    493                                     psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s "
    494                                              "for non-linearity correction --- ignored.\n", concept);
    495                                 } else if (conceptValueItem->type != PS_META_STR) {
    496                                     psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as"
    497                                              " expected for non-linearity correction --- ignored.\n",
    498                                              concept);
     535                            {
     536                                // This is a menu
     537                                psMetadata *options = dataItem->data.V; // Options with concept values as keys
     538                                bool mdok = false; // Success of MD lookup
     539                                psString concept = psMetadataLookupString(&mdok, recipe, "NONLIN.SOURCE");
     540                                if (! mdok || ! concept) {
     541                                    psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction desired, but "
     542                                             "unable to find NONLIN.SOURCE in recipe --- ignored.\n");
    499543                                } else {
    500                                     psString conceptValue = conceptValueItem->data.V;
    501                                     // Get the value of the concept
    502                                     psMetadataItem *optionItem = psMetadataLookup(options, conceptValue);
    503                                     if (!optionItem) {
    504                                         psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA"
    505                                                  " --- ignored.\n", conceptValue);
    506                                     } else if (optionItem->type == PS_META_VECTOR) {
    507                                         // These are the polynomial coefficients
    508                                         psVector *coeff = optionItem->data.V; // The coefficient vector
    509                                         if (coeff->type.type != PS_TYPE_F64) {
    510                                             psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64);
     544                                    psMetadataItem *conceptValueItem = pmCellGetConcept(inputCell, concept);
     545                                    if (! conceptValueItem) {
     546                                        psLogMsg("phase2", PS_LOG_WARN, "Unable to find value of concept %s "
     547                                                 "for non-linearity correction --- ignored.\n", concept);
     548                                    } else if (conceptValueItem->type != PS_META_STR) {
     549                                        psLogMsg("phase2", PS_LOG_WARN, "Type for concept %s isn't STRING, as"
     550                                                 " expected for non-linearity correction --- ignored.\n",
     551                                                 concept);
     552                                    } else {
     553                                        psString conceptValue = conceptValueItem->data.V;
     554                                        // Get the value of the concept
     555                                        psMetadataItem *optionItem = psMetadataLookup(options, conceptValue);
     556                                        if (!optionItem) {
     557                                            psLogMsg("phase2", PS_LOG_WARN, "Unable to find %s in NONLIN.DATA"
     558                                                     " --- ignored.\n", conceptValue);
     559                                        } else if (optionItem->type == PS_META_VEC) {
     560                                            // These are the polynomial coefficients
     561                                            psVector *coeff = optionItem->data.V; // The coefficient vector
     562                                            if (coeff->type.type != PS_TYPE_F64) {
     563                                                psVector *temp = psVectorCopy(NULL, coeff, PS_TYPE_F64);
     564                                                psFree(coeff);
     565                                                coeff = temp;
     566                                            }
     567                                            // Polynomial correction
     568                                            psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,
     569                                                                                             coeff->n + 1);
     570                                            for (int i = 0; i < coeff->n; i++) {
     571                                                correction->coeff[i] = coeff->data.F64[i];
     572                                            }
     573                                            (void)pmNonLinearityPolynomial(inputReadout, correction);
    511574                                            psFree(coeff);
    512                                             coeff = temp;
     575                                            psFree(correction);
     576                                        } else if (optionItem->type == PS_META_STR) {
     577                                            // This is a filename
     578                                            psString tableName = optionItem->data.V; // The filename
     579                                            psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0);
     580                                            int numLines = 0; // Number of lines read from table
     581                                            if ((numLines = psLookupTableRead(table)) <= 0) {
     582                                                psErrorStackPrint(stderr, "Unable to read non-linearity "
     583                                                                  "correction file %s --- ignored\n",
     584                                                                  tableName);
     585                                            } else {
     586#ifdef PRODUCTION
     587                                                (void)pmNonLinearityLookup(inputReadout, table);
     588#else
     589                                                printf("XXX: Non-linearity correction from lookup table not "
     590                                                       "yet implemented.\n");
     591#endif
     592                                            }
     593                                            psFree(table);
     594                                        } else {
     595                                            psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction "
     596                                                     "desired but unable to interpret NONLIN.DATA for %s"
     597                                                     " --- ignored\n", conceptValue);
    513598                                        }
    514                                         // Polynomial correction
    515                                         psPolynomial1D *correction = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,
    516                                                                                          coeff->n + 1);
    517                                         for (int i = 0; i < coeff->n; i++) {
    518                                             correction->coeff[i] = coeff->data.F64[i];
    519                                         }
    520                                         (void)pmNonLinearityPolynomial(inputReadout, correction);
    521                                         psFree(coeffCopy);
    522                                         psFree(coeff);
    523                                         psFree(correction);
    524                                     } else if (optionItem->type == PS_META_STR) {
    525                                         // This is a filename
    526                                         psString tableName = optionItem->data.V; // The filename
    527                                         psLookupTable *table = psLookupTableAlloc(tableName, "%f %f", 0);
    528                                         if ((int numLines = psLookupTableRead(table)) <= 0) {
    529                                             psErrorStackPrint(stderr, "Unable to read non-linearity "
    530                                                               "correction file %s --- ignored\n", tableName);
    531                                         } else {
    532                                             (void)pmNonLinearityLookup(inputReadout, table);
    533                                         }
    534                                         psFree(table);
    535                                     } else {
    536                                         psLogMsg("phase2", PS_LOG_WARN, "Non-linearity correction "
    537                                                  "desired but unable to interpret NONLIN.DATA for %s"
    538                                                  " --- ignored\n", conceptValue);
    539599                                    }
    540600                                }
    … …  
    553613                // it's not in there yet.  Once it is, we can get rid of "pedestal".
    554614
     615#ifdef PRODUCTION
    555616                psImage *pedestal = NULL;       // Pedestal image (bias + scaled dark)
    556617                if (doDark) {
    557                     float inputTime = pmReadoutGetDarkTime(inputReadout); // Dark time for input image
     618                    // Dark time for input image
     619                    float inputTime = psMetadataLookupF32(NULL, inputCell->concepts, "CELL.DARKTIME");
    558620                    if (inputTime <= 0) {
    559621                        psErrorStackPrint(stderr, "DARKTIME for input image (%f) is non-positive.\n",
    … …  
    562624                    }
    563625
    564                     pedestal = psBinaryOp(NULL, darkImage, "*",
    565                                           psScalarAlloc(inputTime * darkTime, PS_TYPE_F32));
     626                    pedestal = (psImage*)psBinaryOp(NULL, darkImage, "*",
     627                                                    psScalarAlloc(inputTime * darkTime, PS_TYPE_F32));
    566628                }
    567629                if (doBias) {
    … …  
    569631                        if (biasImage->numRows != darkImage->numRows ||
    570632                            biasImage->numCols != darkImage->numCols) {
    571                             psError("phase2", true, "Bias and dark images have different dimensions: "
     633                            psError(PS_ERR_IO, true, "Bias and dark images have different dimensions: "
    572634                                    "%dx%d vs %dx%d\n", biasImage->numCols, biasImage->numRows,
    573635                                    darkImage->numCols, darkImage->numRows);
    … …  
    579641                    }
    580642                }
    581 
     643#endif
    582644                if (overscanMode == PM_OVERSCAN_ROWS || overscanMode == PM_OVERSCAN_COLUMNS) {
    583645                    // Need to get the read direction
    584                     int readdir = pmCellGetReaddir(inputCell); // Read direction
     646                    int readdir = psMetadataLookupS32(NULL, inputCell->concepts, "CELL.READDIR");
    585647                    if (readdir == 1) {
    586648                        overscanMode = PM_OVERSCAN_ROWS;
    … …  
    594656                }
    595657
    596                 void *overscanResult = NULL; // Result of overscan fit
    597                 (void)pmSubtractBias(inputImage, overscanResult, inputBias, overscanMode, overscanStat,
    598                                      overscanBin, overscanFit, pedestal);
    599 
    600                 // Output overscan fit results, if required
    601                 if (doOverscan) {
    602                     switch (overscanFit) {
    603                       case PM_FIT_POLYNOMIAL:
    604                         psPolynomial1D *poly = (psPolynomial1D*)overscanResult; // The polynomial
    605                         psString coeffs = NULL; // String containing the coefficients
    606                         for (int i = 0; i < poly->n; i++) {
    607                             psStringAppend(&coeffs, "%.2f ", poly->coeffs[i]);
     658                printf("mode: %d\n", overscanMode);
     659
     660                if (doBias || doOverscan || doDark) {
     661                    void *overscanResult = NULL; // Result of overscan fit
     662#ifdef PRODUCTION
     663                    (void)pmSubtractBias(inputReadout, overscanResult, inputOverscans, overscanMode,
     664                                         overscanStats, overscanBins, overscanFit, pedestal);
     665#else
     666                    (void)pmSubtractBias(inputReadout, overscanResult, inputOverscans, overscanMode,
     667                                         overscanStats, overscanBins, overscanFit, biasReadout);
     668#endif
     669                   
     670                    // Output overscan fit results, if required
     671                    if (doOverscan) {
     672                        switch (overscanFit) {
     673                          case PM_FIT_POLYNOMIAL:
     674                            {
     675                                psPolynomial1D *poly = (psPolynomial1D*)overscanResult; // The polynomial
     676                                psString coeffs = NULL; // String containing the coefficients
     677                                for (int i = 0; i < poly->n; i++) {
     678                                    psStringAppend(&coeffs, "%.2f ", poly->coeff[i]);
     679                                }
     680                                psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n",
     681                                         coeffs);
     682                                psFree(coeffs);
     683                            }
     684                            break;
     685                          case PM_FIT_SPLINE:
     686                            {
     687                                psSpline1D *spline = (psSpline1D*)overscanResult; // The spline
     688                                psString coeffs = NULL; // String containing the coefficients
     689                                for (int i = 0; i < spline->n; i++) {
     690                                    psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
     691                                    psStringAppend(&coeffs, "%d: ", i);
     692                                    for (int j = 0; j < poly->n; j++) {
     693                                        psStringAppend(&coeffs, "%.2f ", poly->coeff[i]);
     694                                    }
     695                                    psStringAppend(&coeffs, "\n");
     696                                }
     697                                psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n",
     698                                         coeffs);
     699                                psFree(coeffs);
     700                            }
     701                            break;
     702                          case PM_FIT_NONE:
     703                            break;
     704                          default:
     705                            psAbort(__func__, "Should never get here!!!\n");
    608706                        }
    609                         psLogMsg("phase2", PS_LOG_INFO, "Overscan polynomial coefficients:\n%s\n",
    610                                  coeffs);
    611                         psFree(coeffs);
    612                         break;
    613                       case PM_FIT_SPLINE:
    614                         psSpline1D *spline = (psSpline1D*)overscanResult; // The spline
    615                         psString coeffs = NULL; // String containing the coefficients
    616                         for (int i = 0; i < spline->n; i++) {
    617                             psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
    618                             psStringAppend(&coeffs, "%d: ", i);
    619                             for (int j = 0; j < poly->n; j++) {
    620                                 psStringAppend(&coeffs, "%.2f ", poly->coeffs[i]);
    621                             }
    622                             psStringAppend(&coeffs, "\n");
    623                         }
    624                         psLogMsg("phase2", PS_LOG_INFO, "Overscan spline coefficients:\n%s\n", coeffs);
    625                         psFree(coeffs);
    626                         break;
    627                       case PM_FIT_NONE:
    628                         break;
    629                       default:
    630                         psAbort("Should never get here!!!\n");
    631707                    }
    632708                }
    … …  
    634710                // Flat-field correction
    635711                if (doFlat) {
    636                     (void)pmFlatField(inputReadout, maskReadout, flatReadout);
     712                    (void)pmFlatField(inputReadout, flatReadout);
    637713                }
    638714
    … …  
    659735
    660736    // Write the output
    661     pmFPAWrite(outputFile, input);
    662     pmFPAWriteMask(outputMaskFile, input);
     737    pmFPAWrite(outputFile, input, database);
     738//    pmFPAWriteMask(outputMaskFile, input);
    663739#ifdef PRODUCTION
    664740    psFitsClose(outputFile);
    665     psFitsClose(outputMaskFile);
     741//    psFitsClose(outputMaskFile);
    666742#else
    667743    psFree(outputFile);
    668     psFree(outputMaskFile);
    669 #endif
    670 
    671 #endif
    672 
     744//    psFree(outputMaskFile);
     745#endif
     746
     747    psFree(arguments);
    673748    psFree(site);
    674749    psFree(header);
    … …  
    681756    psFree(flat);
    682757    psFree(overscanStats);
    683     psFree(inputFile);
    684758
    685759#if 1
Note: See TracChangeset for help on using the changeset viewer.