IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7260 for trunk/ppMerge/src


Ignore:
Timestamp:
Jun 1, 2006, 10:16:15 AM (20 years ago)
Author:
Paul Price
Message:

Working version for flat-fielding. Need to add extenal scale/zero input, mask saturated pixels, test on dark, bias, fringe.

Location:
trunk/ppMerge/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppMerge/src/Makefile.am

    r7068 r7260  
    11bin_PROGRAMS = ppMerge
    22
    3 ppMerge_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(ppMerge_CFLAGS)
    4 ppMerge_LDADD = $(PSLIB_LIBS) $(PSMODULE_LIBS)
     3ppMerge_CFLAGS += -g $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     4ppMerge_LDFLAGS = $(PSMODULE_LIBS) $(PSLIB_LIBS)
     5### For profiling:
     6#ppMerge_CFLAGS += -g -pg -fprofile-arcs $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     7#ppMerge_LDFLAGS = -pg -Wl,--start-group -Wl,-Bstatic $(PSMODULE_LIBS) $(PSLIB_LIBS) -Wl,-Bdynamic
     8
    59ppMerge_SOURCES =               \
    610        ppMerge.c               \
    … …  
    2226        ppMergeOptions.h
    2327
     28CLEANFILES = *~
     29
    2430clean-local:
    2531        -rm -f TAGS

  • trunk/ppMerge/src/ppMerge.c

    r7067 r7260  
    2121int main(int argc, char **argv)
    2222{
     23    psMemThreadSafety(false);           // Turn off thread safety, for more
    2324    psTimerStart(TIMERNAME);
    2425
    … …  
    4748    ppMergeCombine(scale, zero, data, options, config);
    4849
     50    pmFPAPrint(data->out, true, true);
     51
    4952    // Cleaning up
    5053    psFree(data);

  • trunk/ppMerge/src/ppMergeBackground.c

    r7067 r7260  
    4949    // Sanity checks
    5050    assert(!options->scale || scales);
    51     assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F64 &&
     51    assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F32 &&
    5252                                   (*scales)->numCols == data->numCells &&
    5353                                   (*scales)->numRows == filenames->n));
    5454    assert(!options->zero || zeros);
    55     assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F64 &&
     55    assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F32 &&
    5656                                 (*zeros)->numCols == data->numCells &&
    5757                                 (*zeros)->numRows == filenames->n));
    5858
    59     psImage *background = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64); // Background measurements
     59    psImage *background = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); // Background measurements
    6060    psImageInit(background, NAN);
    6161    psStats *bgStats = psStatsAlloc(options->background); // Statistic to measure the background
    … …  
    6363    psImage *exptime = NULL;            // The exposure time for each integration of each cell
    6464    if (options->scale) {
    65         gains = psVectorAlloc(data->numCells, PS_TYPE_F64);
     65        gains = psVectorAlloc(data->numCells, PS_TYPE_F32);
     66        gains->n = data->numCells;
    6667    }
    6768    if (options->exptime) {
    68         exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64);
     69        exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F32);
    6970    }
    7071
    … …  
    7576            continue;
    7677        }
     78        psTrace(__func__, 9, "Opening %s to get background...\n", name);
    7779        psFits *inFile = psFitsOpen(filenames->data[i], "r"); // The FITS file to read
    7880        if (!inFile) {
    … …  
    107109                if (options->exptime) {
    108110                    bool mdok = true;   // Status of MD lookup
    109                     exptime->data.F64[i][cellNum] = psMetadataLookupF32(&mdok, cell->concepts,
     111                    exptime->data.F32[i][cellNum] = psMetadataLookupF32(&mdok, cell->concepts,
    110112                                                                        "CELL.EXPOSURE");
    111                     if (!mdok || isnan(exptime->data.F64[i][cellNum])) {
     113                    if (!mdok || isnan(exptime->data.F32[i][cellNum])) {
    112114                        psLogMsg(__func__, PS_LOG_WARN, "CELL.EXPOSURE for file %s chip %d cell %d is not "
    113115                                 "set.\n", name, j, k);
    114                         exptime->data.F64[i][cellNum] = NAN;
     116                        exptime->data.F32[i][cellNum] = NAN;
    115117                        status = false;
    116118                    }
    … …  
    151153
    152154                    // Get the background
    153                     psImageStats(bgStats, image, readout->mask, options->combine->maskVal);
    154                     background->data.F64[i][cellNum] = getStat(bgStats, options->background);
     155                    int sampleSize = (image->numCols * image->numRows) / options->sample; // Size of sample
     156                    psVector *sample = psVectorAlloc(sampleSize, PS_TYPE_F32); // Sample of the image
     157                    sample->n = sampleSize;
     158                    psVector *sampleMask = NULL; // Mask for sample
     159                    if (readout->mask) {
     160                        sampleMask = psVectorAlloc(sampleSize, PS_TYPE_U8);
     161                        sampleMask->n = sampleSize;
     162                    }
     163                    psImage *mask = readout->mask; // The mask image
     164                    for (long i = 0; i < sampleSize; i++) {
     165                        int j = i * options->sample; // Index into image
     166                        int x = j % image->numCols; // x index
     167                        int y = j / image->numCols; // y index
     168                        sample->data.F32[i] = image->data.F32[y][x];
     169                        if (readout->mask) {
     170                            sampleMask->data.U8[i] = mask->data.U8[y][x];
     171                        }
     172                    }
     173                    psVectorStats(bgStats, sample, sampleMask, NULL, options->combine->maskVal);
     174                    psFree(sample);
     175                    psFree(sampleMask);
     176                    background->data.F32[i][cellNum] = getStat(bgStats, options->background);
     177                    psTrace(__func__, 3, "Background for %s, cell %d is %f\n", name, cellNum,
     178                            background->data.F32[i][cellNum]);
    155179                }
    156180
    … …  
    167191        // Need to normalize over the focal plane
    168192        bool converge = true;           // Did we converge?
     193        if (psTraceGetLevel(__func__) > 9) {
     194            for (int i = 0; i < gains->n; i++) {
     195                psTrace(__func__, 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]);
     196            }
     197        }
    169198        psVector *fluxes = pmFlatNormalize(&converge, gains, background, MAX_ITER, TOLERANCE);
    170199        if (!converge || !fluxes) {
    … …  
    174203        psFree(gains);
    175204
    176         *scales = psImageCopy(*scales, background, PS_TYPE_F64);
     205        *scales = psImageCopy(*scales, background, PS_TYPE_F32);
    177206        psImage *scalesDeref = *scales; // Dereference the pointer
    178207
    179208        for (int i = 0; i < scalesDeref->numRows; i++) {
    180             psF64 bg = fluxes->data.F64[i];
     209            psF32 bg = fluxes->data.F32[i];
    181210            for (int j = 0; j < scalesDeref->numCols; j++) {
    182                 scalesDeref->data.F64[i][j] = bg;
     211                scalesDeref->data.F32[i][j] = bg;
    183212            }
    184213        }
    … …  
    188217        if (!options->scale) {
    189218            // Copy over the exposure times
    190             psImageCopy(*scales, exptime, PS_TYPE_F64);
     219            psImageCopy(*scales, exptime, PS_TYPE_F32);
    191220        } else {
    192221            psBinaryOp(*scales, *scales, "*", exptime);
    … …  
    199228            *zeros = psMemIncrRefCounter(background);
    200229        } else {
    201             *zeros = psImageCopy(*zeros, background, PS_TYPE_F64);
     230            *zeros = psImageCopy(*zeros, background, PS_TYPE_F32);
     231        }
     232    }
     233
     234    // Diagnostic stuff
     235    if (scales && *scales && psTraceGetLevel(__func__) > 9) {
     236        psImage *scalesDeref = *scales; // Dereference the pointer
     237        for (int i = 0; i < scalesDeref->numRows; i++) {
     238            for (int j = 0; j < scalesDeref->numCols; j++) {
     239                psTrace(__func__, 9, "Scale for exposure %d, cell %d is: %f\n", i, j,
     240                        scalesDeref->data.F32[i][j]);
     241            }
     242        }
     243    }
     244    if (zeros && *zeros && psTraceGetLevel(__func__) > 9) {
     245        psImage *zerosDeref = *zeros; // Dereference the pointer
     246        for (int i = 0; i < zerosDeref->numRows; i++) {
     247            for (int j = 0; j < zerosDeref->numCols; j++) {
     248                psTrace(__func__, 9, "Zero for exposure %d, cell %d is: %f\n", i, j,
     249                        zerosDeref->data.F32[i][j]);
     250            }
    202251        }
    203252    }

  • trunk/ppMerge/src/ppMergeCheckInputs.c

    r7073 r7260  
    2727    psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
    2828    assert(filenames);
     29    if (!data->in) {
     30        data->in = psArrayAlloc(filenames->n);
     31        data->in->n = filenames->n;
     32    }
    2933    int numGood = 0;                    // Number of good files
    3034    for (int i = 0; i < filenames->n; i++) {
    … …  
    3337            continue;
    3438        }
     39        psTrace(__func__, 1, "Checking input file %s....\n", name);
    3540        psFits *inFile = psFitsOpen(filenames->data[i], "r"); // The FITS file to read
    3641        if (!inFile) {
    … …  
    4247        }
    4348        psMetadata *header = psFitsReadHeader(NULL, inFile); // The FITS (primary) header
     49        if (psTraceGetLevel(__func__) > 9) {
     50            psTrace(__func__, 9, "Primary header:\n");
     51            psMetadataPrint(header, 9);
     52        }
    4453        psFitsClose(inFile);
    4554
    … …  
    7483        // Use the first valid input as the basis for the output --- including the header
    7584        if (!data->out) {
     85            psTrace(__func__, 5, "Constructing output using %s as a template.\n", name);
    7686            data->out = pmFPAConstruct(config->camera);
    7787            pmFPAview *view = pmFPAAddSourceFromHeader(data->out, header, options->format);
    7888            psFree(view);
     89#if 0
     90            pmFPACopyStructure(data->out, data->in->data[i], 1, 1);
     91#endif
    7992        }
    8093
     94        psTrace(__func__, 3, "%s checks out.\n", name);
    8195        numGood++;
    8296    }
    … …  
    99113    }
    100114    data->numCells = numCells;
     115    psTrace(__func__, 3, "Output has %d cells.\n", numCells);
    101116
     117    psTrace(__func__, 3, "We have %d good inputs.\n", numGood);
    102118    if (numGood > 1) {
    103119        return data;

  • trunk/ppMerge/src/ppMergeCombine.c

    r7067 r7260  
    77#include "ppMergeData.h"
    88#include "ppMergeCombine.h"
     9
    910
    1011// Combine the inputs
    … …  
    2223    // Sanity checks
    2324    assert(!options->scale || scales);
    24     assert(!scales || (scales->type.type == PS_TYPE_F64 &&
     25    assert(!scales || (scales->type.type == PS_TYPE_F32 &&
    2526                       scales->numCols == data->numCells &&
    2627                       scales->numRows == filenames->n));
    2728    assert(!options->zero || zeros);
    28     assert(!zeros || (zeros->type.type == PS_TYPE_F64 &&
     29    assert(!zeros || (zeros->type.type == PS_TYPE_F32 &&
    2930                      zeros->numCols == data->numCells &&
    3031                      zeros->numRows == filenames->n));
    … …  
    3233    // Iterate over the FPA
    3334    int cellNum = -1;                   // Cell number
    34     pmFPAWrite(data->out, data->outFile, config->database, false);
     35    pmFPAWrite(data->out, data->outFile, config->database, false, false); // Write header only
    3536    psArray *chips = data->out->chips;  // Array of output chips
    3637    for (int i = 0; i < chips->n; i++) {
    … …  
    3940            continue;
    4041        }
    41         pmChipWrite(chip, data->outFile, config->database, false);
     42        pmChipWrite(chip, data->outFile, config->database, false, false); // Write header only
    4243        psArray *cells = chip->cells;   // Array of output cells
    4344        for (int j = 0; j < cells->n; j++) {
    … …  
    4748            }
    4849            cellNum++;
    49             pmCellWrite(cell, data->outFile, config->database, false);
     50            pmCellWrite(cell, data->outFile, config->database, false); // Write header only
    5051            pmReadout *readout = pmReadoutAlloc(cell); // Output readout of interest
    5152            psArray *stack = psArrayAlloc(filenames->n); // Stack of readouts to combine
     53            stack->n = filenames->n;
    5254            psVector *cellScales = NULL; // Scales for this cell
    5355            if (scales) {
    … …  
    5961            }
    6062
    61             bool stillReadingRows = false;  // Still reading rows from any of the files?
    62 
     63            int numRead;  // Number of inputs read
    6364            do {
     65                numRead = 0;
    6466                for (int k = 0; k < filenames->n; k++) {
    6567                    if (! filenames->data[k] || strlen(filenames->data[k]) == 0) {
    … …  
    8284
    8385                    // Only reading and writing the first readout in each cell (plane 0)
    84                     stillReadingRows |= pmReadoutReadNext(stack->data[k], fits, 0, options->rows);
     86                    if (pmReadoutReadNext(stack->data[k], fits, 0, options->rows)) {
     87                        numRead++;
     88                    }
    8589                    psFitsClose(fits);
    8690                }
    8791
    88                 pmReadoutCombine(readout, stack, cellZeros, cellScales, options->combine);
    89                 pmReadoutWriteNext(readout, data->outFile, 0);
     92                if (numRead > 0) {
     93                    pmReadoutCombine(readout, stack, cellZeros, cellScales, options->combine);
     94                    psTrace(__func__, 5, "Chip %d, cell %d, readout size %dx%d\n", i, j,
     95                            readout->image->numCols, readout->image->numRows);
     96                }
    9097
    91             } while (stillReadingRows);
    92 
    93             // Write the pixels
    94             pmCellWrite(cell, data->outFile, config->database, true);
     98            } while (numRead > 0);
    9599
    96100            // Blow away the cell data
    … …  
    101105                pmCellFreeData(cellIn);
    102106            }
    103             pmCellFreeData(cell);
     107
     108            // Write the pixels
     109            if (cell->hdu) {
     110                psTrace(__func__, 5, "Writing out cell HDU.\n");
     111                pmCellWrite(cell, data->outFile, config->database, true);
     112                pmCellFreeData(cell);
     113            }
    104114        }
    105 
    106         // Write the pixels
    107         pmChipWrite(chip, data->outFile, config->database, true);
    108115
    109116        // Blow away the chip data
    … …  
    113120            pmChipFreeData(chipIn);
    114121        }
    115         pmChipFreeData(chip);
    116122
     123        // Write the pixels
     124        if (chip->hdu) {
     125            psTrace(__func__, 5, "Writing out chip HDU.\n");
     126            pmChipWrite(chip, data->outFile, config->database, true, false);
     127            pmChipFreeData(chip);
     128        }
    117129    }
    118130
    119     // Write the pixels
    120     pmFPAWrite(data->out, data->outFile, config->database, true);
     131    if (data->out->hdu) {
     132        // Write the pixels
     133        psTrace(__func__, 5, "Writing out FPA HDU.\n");
     134        pmFPAWrite(data->out, data->outFile, config->database, true, false);
     135    }
    121136
    122137    // Blow away the FPA data

  • trunk/ppMerge/src/ppMergeOptions.c

    r7073 r7260  
    137137    OPTION_PARSE(options->combine->fracLow,  recipe, "FRACLOW",   F32 );
    138138    OPTION_PARSE(options->combine->nKeep,    recipe, "NKEEP",     S32 );
    139     OPTION_PARSE(options->combine->maskVal,  recipe, "MASKVAL",   U8 );
     139    OPTION_PARSE(options->combine->maskVal,  recipe, "MASKVAL",   S32 );
    140140    options->combine->combine = parseStat(recipe, "COMBINE");
    141141    options->background       = parseStat(recipe, "BACKGROUND");
    … …  
    145145    // Set options based on the type of calibration frame
    146146    const char *type = psMetadataLookupStr(NULL, config->arguments, "-type"); // The type of calibration frame
    147     if (strcasecmp(type, "BIAS") == 0) {
     147    if (strlen(type) > 0) {
     148        if (strcasecmp(type, "BIAS") == 0) {
     149            options->zero = false;
     150            options->scale = false;
     151            options->exptime = false;
     152        } else if (strcasecmp(type, "DARK") == 0) {
     153            options->zero = false;
     154            options->scale = false;
     155            options->exptime = true;
     156        } else if (strcasecmp(type, "FLAT") == 0) {
     157            options->zero = false;
     158            options->scale = true;
     159            options->exptime = false;
     160        } else if (strcasecmp(type, "FRINGE") == 0) {
     161            options->zero = true;
     162            options->scale = true;
     163            options->exptime = false;
     164        } else {
     165            psLogMsg(__func__, PS_LOG_WARN, "Unrecognised image type: %s --- assuming BIAS.\n", type);
     166            options->zero = false;
     167            options->scale = false;
     168            options->exptime = false;
     169        }
     170    } else {
     171        psLogMsg(__func__, PS_LOG_WARN, "No calibration type specified; assuming BIAS.\n");
    148172        options->zero = false;
    149173        options->scale = false;
    150174        options->exptime = false;
    151     } else if (strcasecmp(type, "DARK") == 0) {
    152         options->zero = false;
    153         options->scale = false;
    154         options->exptime = true;
    155     } else if (strcasecmp(type, "FLAT") == 0) {
    156         options->zero = false;
    157         options->scale = true;
    158         options->exptime = false;
    159     } else if (strcasecmp(type, "FRINGE") == 0) {
    160         options->zero = true;
    161         options->scale = true;
    162         options->exptime = false;
    163     } else {
    164         psLogMsg(__func__, PS_LOG_WARN, "Unrecognised image type: %s --- ignored.\n", type);
    165175    }
    166176
     177
     178#if 0
    167179    // Or you can set them individually
    168180    OPTION_PARSE(options->zero,    config->arguments, "-zero",    Bool);
    169181    OPTION_PARSE(options->scale,   config->arguments, "-scale",   Bool);
    170182    OPTION_PARSE(options->exptime, config->arguments, "-exptime", Bool);
     183#endif
    171184
    172185    // Number of on/off images
Note: See TracChangeset for help on using the changeset viewer.