IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27427


Ignore:
Timestamp:
Mar 23, 2010, 9:02:41 PM (16 years ago)
Author:
Paul Price
Message:

Something was missing from the merge (probably because I tool r27400 out), so had to get it back.

Location:
trunk
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r27402 r27427  
    6969                               const psVector *mask, // Mask for input readouts
    7070                               const psVector *weightings, // Weighting factors for each image
     71                               const psVector *exposures,  // Exposure time for each image
    7172                               const psVector *addVariance // Additional variance for rejection
    7273    );
    … …  
    8384bool ppStackReadoutFinal(const pmConfig *config,   // Configuration
    8485                         pmReadout *outRO,   // Output readout
     86                         pmReadout *expRO,   // Exposure readout
    8587                         const psArray *readouts, // Input readouts
    8688                         const psVector *mask, // Mask for input readouts
    8789                         const psArray *rejected, // Array with pixels rejected in each image
    8890                         const psVector *weightings, // Weighting factors for each image
     91                         const psVector *exposures,  // Exposure times for each image
    8992                         const psVector *addVariance, // Additional variance for rejection
    9093                         bool safety,                 // Enable safety switch?

  • trunk/ppStack/src/ppStackCleanup.c

    r27402 r27427  
    8686        options->outRO = NULL;
    8787
     88        options->expRO->data_exists = false;
     89        options->expRO->parent->data_exists = false;
     90        options->expRO->parent->parent->data_exists = false;
     91        psFree(options->expRO);
     92        options->expRO = NULL;
     93
    8894        pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
    8995        view->chip = view->cell = 0;        // pmFPAviewFreeData doesn't want to deal with readouts

  • trunk/ppStack/src/ppStackCombineFinal.c

    r27403 r27427  
    1414//#define TESTING                         // Enable test output
    1515
    16 bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances,
    17                          ppStackOptions *options, pmConfig *config, bool safe, bool normalise, bool grow)
     16bool ppStackCombineFinal(ppStackThreadData *stack, psArray *covariances, ppStackOptions *options,
     17                         pmConfig *config, bool safe, bool normalise, bool grow)
    1818{
    1919    psAssert(stack, "Require stack");
    … …  
    2121    psAssert(config, "Require configuration");
    2222
    23     int numCols = target->image->numCols, numRows = target->image->numRows; // Size of image
     23    pmReadout *outRO = options->outRO;                                      // Output readout
     24    pmReadout *expRO = options->expRO;                                      // Exposure readout
     25    int numCols = outRO->image->numCols, numRows = outRO->image->numRows; // Size of image
    2426
    2527    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    … …  
    4648    }
    4749
    48     if (!target->mask) {
    49         target->mask = psImageAlloc(target->image->numCols, target->image->numRows, PS_TYPE_IMAGE_MASK);
     50    if (!outRO->mask) {
     51        outRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     52    }
     53    if (!expRO->mask) {
     54        expRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    5055    }
    5156
    … …  
    6570        }
    6671
    67         // call: ppStackReadoutFinal(config, target, readouts, rejected)
     72        // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
    6873        psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
    69         psArrayAdd(job->args, 1, target);
    7074        psArrayAdd(job->args, 1, thread);
    7175        psArrayAdd(job->args, 1, reject);
    … …  
    108112        }
    109113        if (sumWeights > 0.0) {
    110             target->covariance = psImageCovarianceSum(covariances);
    111             psBinaryOp(target->covariance->image, target->covariance->image, "/",
     114            outRO->covariance = psImageCovarianceSum(covariances);
     115            psBinaryOp(outRO->covariance->image, outRO->covariance->image, "/",
    112116                       psScalarAlloc(sumWeights, PS_TYPE_F32));
    113             psImageCovarianceTransfer(target->variance, target->covariance);
     117            psImageCovarianceTransfer(outRO->variance, outRO->covariance);
    114118        }
    115119    } else {
    116         target->covariance = psImageCovarianceNone();
     120        outRO->covariance = psImageCovarianceNone();
    117121    }
    118122
    … …  
    130134            wcsDone = true;
    131135            pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
    132             pmHDU *outHDU = pmHDUFromCell(target->parent); // Output HDU
    133             pmChip *outChip = target->parent->parent; // Output chip
     136            pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
     137            pmChip *outChip = outRO->parent->parent; // Output chip
    134138            pmFPA *outFPA = outChip->parent; // Output FPA
    135139            if (!outHDU || !inHDU) {
    … …  
    154158    }
    155159
     160    // Set exposure time correctly
     161    {
     162        float exptime = 0.0;            // Summed exposure time
     163        for (int i = 0; i < options->num; i++) {
     164            if (options->inputMask) {
     165                continue;
     166            }
     167            exptime += options->exposures->data.F32[i];
     168        }
     169
     170        {
     171            psMetadataItem *item = psMetadataLookup(outRO->parent->concepts, "CELL.EXPOSURE");
     172            item->data.F32 = exptime;
     173        }
     174        {
     175            psMetadataItem *item = psMetadataLookup(outRO->parent->parent->parent->concepts, "FPA.EXPOSURE");
     176            item->data.F32 = exptime;
     177        }
     178    }
     179
    156180    // Put version information into the header
    157     pmHDU *hdu = pmHDUFromCell(target->parent);
     181    pmHDU *hdu = pmHDUFromCell(outRO->parent);
    158182    if (!hdu) {
    159183        psError(PPSTACK_ERR_PROG, false, "Unable to find HDU for output.");
    … …  
    171195    psStringAppend(&name, "combined_image_final_%d.fits", pass);
    172196    pass++;
    173     ppStackWriteImage(name, NULL, target->image, config);
     197    ppStackWriteImage(name, NULL, outRO->image, config);
    174198    psStringSubstitute(&name, "mask", "image");
    175     ppStackWriteImage(name, NULL, target->mask, config);
     199    ppStackWriteImage(name, NULL, outRO->mask, config);
    176200    psStringSubstitute(&name, "variance", "mask");
    177     ppStackWriteImage(name, NULL, target->variance, config);
     201    ppStackWriteImage(name, NULL, outRO->variance, config);
    178202    psFree(name);
    179203
    180     pmStackVisualPlotTestImage(target->image, "combined_image_final.fits");
     204    pmStackVisualPlotTestImage(outRO->image, "combined_image_final.fits");
    181205#endif
    182206

  • trunk/ppStack/src/ppStackCombineInitial.c

    r27402 r27427  
    8989    }
    9090
     91    ppStackMemDump("initial");
     92
    9193#ifdef TESTING
    9294    ppStackWriteImage("combined_image_initial.fits", NULL, options->outRO->image, config);

  • trunk/ppStack/src/ppStackCombinePrepare.c

    r27402 r27427  
    1010#include "ppStackLoop.h"
    1111
    12 bool ppStackCombinePrepare(pmReadout **ro, const char *name, ppStackFileList files, ppStackThreadData *stack,
     12bool ppStackCombinePrepare(const char *outName, const char *expName,
     13                           ppStackFileList files, ppStackThreadData *stack,
    1314                           ppStackOptions *options, pmConfig *config)
    1415{
    … …  
    3233    }
    3334
    34     pmCell *cell = pmFPAfileThisCell(config->files, view, name); // Output cell
    35     *ro = pmReadoutAlloc(cell); // Output readout
     35    pmCell *cell = pmFPAfileThisCell(config->files, view, outName); // Output cell
     36    options->outRO = pmReadoutAlloc(cell); // Output readout
     37
     38    pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell
     39    options->expRO = pmReadoutAlloc(expCell); // Output readout
    3640
    3741    psFree(view);
    … …  
    4246    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    4347
    44     if (!pmReadoutStackDefineOutput(*ro, col0, row0, numCols, numRows, true, true, maskBad)) {
     48    if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
     49        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
     50        return false;
     51    }
     52
     53    if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) {
    4554        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
    4655        return false;

  • trunk/ppStack/src/ppStackLoop.c

    r27402 r27427  
    111111
    112112    // Prepare for combination
    113     if (!ppStackCombinePrepare(&options->outRO, "PPSTACK.OUTPUT", PPSTACK_FILES_STACK,
     113    if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK,
    114114                               stack, options, config)) {
    115115        psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
    … …  
    160160    // Final combination
    161161    psTrace("ppStack", 2, "Final stack of convolved images....\n");
    162     if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config,
    163                              false, false, true)) {
     162    if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true)) {
    164163        psError(psErrorCodeLast(), false, "Unable to perform final combination.");
    165164        psFree(stack);
    … …  
    203202
    204203        // Prepare for combination
    205         if (!ppStackCombinePrepare(&options->unconvRO, "PPSTACK.UNCONV", PPSTACK_FILES_UNCONV,
     204        if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV,
    206205                                   stack, options, config)) {
    207206            psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
    … …  
    211210
    212211        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
    213         if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config,
     212        if (!ppStackCombineFinal(stack, options->origCovars, options, config,
    214213                                 false, true, false)) {
    215214            psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination.");
    … …  
    226225        }
    227226        ppStackFileActivation(config, PPSTACK_FILES_UNCONV, false);
    228         options->unconvRO->data_exists = false;
    229         options->unconvRO->parent->data_exists = false;
    230         options->unconvRO->parent->parent->data_exists = false;
    231         psFree(options->unconvRO);
    232         options->unconvRO = NULL;
     227        options->outRO->data_exists = false;
     228        options->outRO->parent->data_exists = false;
     229        options->outRO->parent->parent->data_exists = false;
     230        psFree(options->outRO);
     231        options->outRO = NULL;
     232        options->expRO->data_exists = false;
     233        options->expRO->parent->data_exists = false;
     234        options->expRO->parent->parent->data_exists = false;
     235        psFree(options->expRO);
     236        options->expRO = NULL;
    233237
    234238        for (int i = 0; i < options->num; i++) {

  • trunk/ppStack/src/ppStackLoop.h

    r27402 r27427  
    3838// Prepare for combination
    3939bool ppStackCombinePrepare(
    40     pmReadout **readout,                // Readout to set
    41     const char *name,                   // Name of file
     40    const char *outName,                // Name of output file
     41    const char *expName,                // Name of exposure file
    4242    ppStackFileList files,              // Files of interest
    4343    ppStackThreadData *stack,           // Stack
    … …  
    6161// Final combination
    6262bool ppStackCombineFinal(
    63     pmReadout *target,                  // Target readout
    6463    ppStackThreadData *stack,           // Stack
    6564    psArray *covariances,               // Covariances

  • trunk/ppStack/src/ppStackOptions.c

    r27402 r27427  
    2121    psFree(options->psf);
    2222    psFree(options->inputSeeing);
     23    psFree(options->exposures);
    2324    psFree(options->inputMask);
    2425    psFree(options->sourceLists);
    … …  
    3233    psFree(options->convCovars);
    3334    psFree(options->outRO);
    34     psFree(options->unconvRO);
     35    psFree(options->expRO);
    3536    psFree(options->inspect);
    3637    psFree(options->rejected);
    … …  
    6162    options->zp = NAN;
    6263    options->inputSeeing = NULL;
     64    options->exposures = NULL;
    6365    options->targetSeeing = NAN;
    6466    options->inputMask = NULL;
    … …  
    7577    options->convCovars = NULL;
    7678    options->outRO = NULL;
    77     options->unconvRO = NULL;
     79    options->expRO = NULL;
    7880    options->inspect = NULL;
    7981    options->rejected = NULL;

  • trunk/ppStack/src/ppStackOptions.h

    r27402 r27427  
    2222    psVector *inputSeeing;              // Input seeing FWHMs
    2323    float targetSeeing;                 // Target seeing FWHM
     24    psVector *exposures;                // Exposure times
    2425    float sumExposure;                  // Sum of exposure times
    2526    float zp;                           // Zero point for output
    … …  
    3839    // Combine initial
    3940    pmReadout *outRO;                   // Output readout
    40     pmReadout *unconvRO;                // Unconvolved readout
     41    pmReadout *expRO;                   // Exposure readout
    4142    psArray *inspect;                   // Array of arrays of pixels to inspect
    4243    // Rejection

  • trunk/ppStack/src/ppStackPrepare.c

    r27402 r27427  
    126126    options->inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
    127127    psVectorInit(options->inputMask, 0);
     128    options->exposures = psVectorAlloc(options->num, PS_TYPE_F32);
     129    psVectorInit(options->exposures, NAN);
    128130
    129131    pmFPAfileActivate(config->files, false, NULL);
    … …  
    134136    }
    135137
    136     psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, "^PPSTACK.INPUT$");
    137     psMetadataItem *fileItem; // Item from iteration
    138138    psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
    139139    int numCols = 0, numRows = 0;   // Size of image
    140     int index = 0;                  // Index for file
    141     while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
    142         assert(fileItem->type == PS_DATA_UNKNOWN);
    143         pmFPAfile *inputFile = fileItem->data.V; // An input file
     140    options->sumExposure = 0.0;
     141    for (int i = 0; i < num; i++) {
     142        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
     143        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
     144
     145        options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE");
     146        options->sumExposure += options->exposures->data.F32[i];
    144147
    145148        // Get list of PSFs, to determine target PSF
    146149        if (options->convolve) {
    147             pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
     150            pmChip *chip = pmFPAviewThisChip(view, file->fpa); // The chip: holds the PSF
    148151            pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
    149152            if (!psf) {
    150153                psError(PPSTACK_ERR_PROG, false, "Unable to find PSF.");
    151154                psFree(view);
    152                 psFree(fileIter);
    153155                psFree(psfs);
    154156                return false;
    155157            }
    156             psfs->data[index] = psMemIncrRefCounter(psf);
     158            psfs->data[i] = psMemIncrRefCounter(psf);
    157159            psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF");
    158160
    159             pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
     161            pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
    160162            pmHDU *hdu = pmHDUFromCell(cell);
    161163            assert(hdu && hdu->header);
    … …  
    165167                psError(PPSTACK_ERR_PROG, false, "Unable to determine size of image from PSF.");
    166168                psFree(view);
    167                 psFree(fileIter);
    168169                psFree(psfs);
    169170                return false;
    … …  
    180181        pmDetections *detections = NULL;
    181182        if (options->convolve || options->matchZPs || options->photometry || redoPhot) {
    182             pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
     183            pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout with sources
    183184            detections = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.DETECTIONS"); // Sources
    184185            if (!detections) {
    … …  
    188189            psAssert (detections->allSources, "missing sources?");
    189190
    190             options->sourceLists->data[index] = psMemIncrRefCounter(detections->allSources);
     191            options->sourceLists->data[i] = psMemIncrRefCounter(detections->allSources);
    191192        }
    192193
    193194        // Re-do photometry if we don't trust the source lists
    194195        if (redoPhot) {
    195             psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
     196            psTrace("ppStack", 2, "Photometering input %d of %d....\n", i, num);
    196197            pmFPAfileActivate(config->files, false, NULL);
    197             ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
     198            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i);
    198199            if (options->convolve) {
    199200                pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
    200201            }
    201             pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
     202            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File
    202203            pmFPAview *photView = ppStackFilesIterateDown(config);
    203204            if (!photView) {
    … …  
    223224            ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
    224225        }
    225 
    226         index++;
    227     }
    228     psFree(fileIter);
     226    }
    229227
    230228    psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message

  • trunk/ppStack/src/ppStackReadout.c

    r27402 r27427  
    2323    psVector *mask = options->inputMask; // Mask for inputs
    2424    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     25    psVector *exposures = options->exposures;   // Exposure times for each image
    2526    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
    2627
    2728    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
    28                                          weightings, addVariance);
     29                                         weightings, exposures, addVariance);
    2930    thread->busy = false;
    3031
    … …  
    3738
    3839    psArray *args = job->args;          // Arguments
    39     pmReadout *target = args->data[0];  // Output readout
    40     ppStackThread *thread = args->data[1]; // Thread
    41     psArray *reject = args->data[2];    // Rejected pixels for each image
    42     ppStackOptions *options = args->data[3]; // Options
    43     pmConfig *config = args->data[4];   // Configuration
    44     bool safety = PS_SCALAR_VALUE(args->data[5], U8);    // Safety switch on?
    45     bool normalise = PS_SCALAR_VALUE(args->data[6], U8); // Normalise images?
     40    ppStackThread *thread = args->data[0]; // Thread
     41    psArray *reject = args->data[1];    // Rejected pixels for each image
     42    ppStackOptions *options = args->data[2]; // Options
     43    pmConfig *config = args->data[3];   // Configuration
     44    bool safety = PS_SCALAR_VALUE(args->data[4], U8);    // Safety switch on?
     45    bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images?
    4646
    4747    psVector *mask = options->inputMask; // Mask for inputs
    4848    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     49    psVector *exposures = options->exposures;   // Exposure times for each image
    4950    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
    5051    psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
    5152
    52     bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, reject,
    53                                       weightings, addVariance, safety, norm); // Status of operation
     53    bool status = ppStackReadoutFinal(config, options->outRO, options->expRO, thread->readouts, mask, reject,
     54                                      weightings, exposures, addVariance, safety, norm); // Status of operation
    5455
    5556    thread->busy = false;
     57
     58    psAssert(status, "Stacking failed.");
    5659
    5760    return status;
    … …  
    101104
    102105psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    103                                const psVector *mask, const psVector *weightings, const psVector *addVariance)
     106                               const psVector *mask, const psVector *weightings, const psVector *exposures,
     107                               const psVector *addVariance)
    104108{
    105109    assert(config);
    … …  
    149153        }
    150154
    151         stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
    152     }
    153 
    154     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
     155        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
     156                                          addVariance->data.F32[i]);
     157    }
     158
     159    if (!pmStackCombine(outRO, NULL, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
    155160                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
    156161        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    … …  
    187192
    188193
    189 bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
     194bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts,
    190195                         const psVector *mask, const psArray *rejected, const psVector *weightings,
    191                          const psVector *addVariance, bool safety, const psVector *norm)
     196                         const psVector *exposures, const psVector *addVariance, bool safety,
     197                         const psVector *norm)
    192198{
    193199    assert(config);
    194200    assert(outRO);
     201    assert(expRO);
    195202    assert(readouts);
    196203    assert(!rejected || readouts->n == rejected->n);
    … …  
    238245        }
    239246
    240         pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i],
     247        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
    241248                                             addVariance ? addVariance->data.F32[i] : NAN);
    242249        data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
    … …  
    250257    }
    251258
    252     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
     259    if (!pmStackCombine(outRO, expRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
    253260                        combineSys, combineDiscard, useVariance, safe, rejected)) {
    254261        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    … …  
    259266    pmCell *outCell = outRO->parent;    // Output cell
    260267    pmChip *outChip = outCell->parent;  // Output chip
    261 
    262268    outRO->data_exists = true;
    263269    outCell->data_exists = true;
    264270    outChip->data_exists = true;
    265271
     272    pmCell *expCell = expRO->parent;    // Exposure cell
     273    pmChip *expChip = expCell->parent;  // Exposure chip
     274    expRO->data_exists = true;
     275    expCell->data_exists = true;
     276    expChip->data_exists = true;
     277
    266278    psFree(stack);
    267279

  • trunk/ppStack/src/ppStackReject.c

    r27426 r27427  
    162162        psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
    163163                 psTimerClear("PPSTACK_REJECT"));
    164 
    165         ppStackMemDump("reject");
    166164    }
    167165

  • trunk/ppStack/src/ppStackSources.c

    r27402 r27427  
    6464
    6565    if (!options->matchZPs && !options->photometry) {
    66         int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
    67         options->norm = psVectorAlloc(num, PS_TYPE_F32);
    68         psVectorInit (options->norm, 0.0);
    69 
    70         // XXX do I need to set this?
    71         // options->sumExposure = sumExpTime;
    72 
     66        options->norm = psVectorAlloc(options->num, PS_TYPE_F32);
     67        psVectorInit(options->norm, 0.0);
    7368        return true;
    7469    }
    … …  
    137132    }
    138133
    139     int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
     134    int num = options->num;            // Number of inputs
    140135    psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
    141136
    … …  
    146141    float airmassTerm = NAN;            // Airmass term
    147142    float zpTarget = NAN;               // Target zero point
    148     float sumExpTime = 0.0;             // Sum of the exposure time
    149143    int numGoodImages = 0;              // Number of good images
    150144    for (int i = 0; i < num; i++) {
    … …  
    160154
    161155        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
    162         pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
    163156
    164157#if defined(TESTING) && 0
    … …  
    177170#endif
    178171
    179         float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
     172        float exptime = options->exposures->data.F32[i]; // Exposure time
    180173        float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass
    181174        const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
    … …  
    221214
    222215        zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
    223         sumExpTime += exptime;
    224 
    225     }
    226 
    227     options->sumExposure = sumExpTime;
     216    }
    228217
    229218    if (numGoodImages == 0) {
    … …  
    291280                }
    292281                psArray *sources = sourceLists->data[i]; // Sources of interest
    293                 float magCorr = zp->data.F32[i] + trans->data.F32[i] - 2.5*log10(sumExpTime);
     282                float magCorr = zp->data.F32[i] + trans->data.F32[i] - 2.5*log10(options->sumExposure);
    294283                if (zpExpNum == numGoodImages) {
    295284                    // Using measured zero points, so attempt to set target zero point

  • trunk/ppStack/src/ppStackThread.c

    r27402 r27427  
    284284
    285285    {
    286         psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7);
     286        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6);
    287287        task->function = &ppStackReadoutFinalThread;
    288288        psThreadTaskAdd(task);

  • trunk/psModules/src/imcombine/pmStack.c

    r27420 r27427  
    3535
    3636//#define TESTING                         // Enable test output
    37 //#define TEST_X 2148-1                     // x coordinate to examine
    38 //#define TEST_Y 248-1                     // y coordinate to examine
     37//#define TEST_X 843-1                     // x coordinate to examine
     38//#define TEST_Y 813-1                     // y coordinate to examine
    3939//#define TEST_RADIUS 0                    // Radius to examine
    4040
    … …  
    4646    psVector *variances;                // Pixel variances
    4747    psVector *weights;                  // Pixel weightings
     48    psVector *exps;                     // Pixel exposures
    4849    psVector *sources;                  // Pixel sources (which image did they come from?)
    4950    psVector *limits;                   // Rejection limits
    … …  
    5758    psFree(buffer->variances);
    5859    psFree(buffer->weights);
     60    psFree(buffer->exps);
    5961    psFree(buffer->sources);
    6062    psFree(buffer->limits);
    … …  
    7375    buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
    7476    buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
     77    buffer->exps = psVectorAlloc(numImages, PS_TYPE_F32);
    7578    buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
    7679    buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
    … …  
    9598static bool combinationMeanVariance(float *mean, // Mean value, to return
    9699                                    float *var, // Variance value, to return
     100                                    float *exp, // Exposure time, to return
     101                                    float *expWeight,          // Weighted exposure time, to return
    97102                                    const psVector *values, // Values to combine
    98103                                    const psVector *variances, // Pixel variances to combine
     104                                    const psVector *exps,      // Exposure times to combine
    99105                                    const psVector *weights // Weights to apply
    100106                                    )
    … …  
    121127    float sumVarianceWeight = 0.0;     // Sum of the pixel variances multiplied by the global weights
    122128    float sumWeight = 0.0;              // Sum of the image weights
     129    float sumExp = 0.0;                 // Sum of the exposure time
     130    float sumExpWeight = 0.0;           // Sum of the exposure time multiplied by the global weights
     131    int numGood = 0;                    // Number of good exposures
    123132    for (int i = 0; i < values->n; i++) {
    124133        sumValueWeight += values->data.F32[i] * weights->data.F32[i];
    … …  
    127136            sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
    128137        }
     138        if (exps) {
     139            sumExp += exps->data.F32[i];
     140            sumExpWeight += exps->data.F32[i] * weights->data.F32[i];
     141            numGood++;
     142        }
    129143    }
    130144
    … …  
    136150    if (var) {
    137151        *var = sumVarianceWeight / PS_SQR(sumWeight);
     152    }
     153    if (exp) {
     154        *exp = sumExp;
     155    }
     156    if (expWeight) {
     157        *expWeight = sumExpWeight;
    138158    }
    139159    return true;
    … …  
    276296                           const psArray *inputs, // Stack data
    277297                           const psVector *weights, // Global (single value) weights for data, or NULL
     298                           const psVector *exps,    // Exposures for data, or NULL
    278299                           const psVector *addVariance, // Additional variance for data
    279300                           const psVector *reject, // Indices of pixels to reject, or NULL
    … …  
    292313    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    293314    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     315    psVector *pixelExps = buffer->exps;       // Exposure times
    294316    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    295317    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
    … …  
    331353        }
    332354        pixelWeights->data.F32[numGood] = data->weight;
     355        pixelExps->data.F32[numGood] = data->exp;
    333356        pixelSources->data.U16[numGood] = i;
    334357        numGood++;
    … …  
    347370    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    348371        for (int i = 0; i < numGood; i++) {
    349             fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %d\n",
     372            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d\n",
    350373                    i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
    351                     addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]);
     374                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
     375                    pixelSuspects->data.U8[i]);
    352376        }
    353377    }
    … …  
    362386                          psImage *mask, // Combined mask, for output
    363387                          psImage *variance, // Combined variance map, for output
     388                          psImage *exp,   // Exposure map (time), for output
     389                          psImage *expnum,       // Exposure map (number) for output
     390                          psImage *expweight,    // Exposure map (weighted time) for output
    364391                          int num,      // Number of good pixels
    365392                          combineBuffer *buffer, // Buffer with vectors
    366393                          int x, int y, // Coordinates of interest; frame of output image
    367394                          psImageMaskType bad, // Value for bad pixels
    368                           bool safe             // Safe combination?
     395                          bool safe,           // Safe combination?
     396                          float invTotalWeight    // Inverse of total weight for all inputs
    369397                          )
    370398{
    … …  
    372400    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    373401    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     402    psVector *pixelExps = buffer->exps;       // Exposure times
    374403
    375404    // Default option is that the pixel is bad
    376405    float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
    377406    psImageMaskType maskValue = bad;    // Value for combined mask
     407    float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted)
     408
    378409    switch (num) {
    379410      case 0: {
    … …  
    393424                  varianceValue = pixelVariances->data.F32[0];
    394425              }
     426              if (exp) {
     427                  expValue = pixelExps->data.F32[0];
     428                  expWeightValue = pixelExps->data.F32[0];
     429              }
    395430              maskValue = 0;
    396431#ifdef TESTING
    … …  
    411446          // Automatically accept the mean of the pixels only if we're not playing safe
    412447          if (!safe) {
    413               float mean, var;   // Mean and variance from combination
    414               if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
    415                   imageValue = mean;
    416                   varianceValue = var;
     448              if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
     449                                          pixelData, pixelVariances, pixelExps, pixelWeights)) {
    417450                  maskValue = 0;
    418451#ifdef TESTING
    419452                  if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    420453                      fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n",
    421                               x, y, mean, var);
     454                              x, y, imageValue, varianceValue);
    422455                  }
    423456#endif
    … …  
    435468      default: {
    436469          // Can combine without too much worrying
    437           float mean, var;           // Mean and variance of the combination
    438           if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
     470          if (!combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
     471                                       pixelData, pixelVariances, pixelExps, pixelWeights)) {
    439472              break;
    440473          }
    441           imageValue = mean;
    442           varianceValue = var;
    443474          maskValue = 0;
    444475#ifdef TESTING
    445476          if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    446               fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, mean, var);
     477              fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
    447478          }
    448479#endif
    … …  
    456487        variance->data.F32[y][x] = varianceValue;
    457488    }
    458 
    459 #ifdef TESTING
    460     mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
    461 #endif
    462 
     489    if (exp) {
     490        exp->data.F32[y][x] = expValue;
     491    }
     492    if (expnum) {
     493        expnum->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
     494    }
     495    if (expweight) {
     496        expweight->data.F32[y][x] = expWeightValue * invTotalWeight;
     497    }
    463498
    464499    return;
    … …  
    803838                              int *numCols, int *numRows, // Size of (sky) images
    804839                              const psArray *input, // Input array of pmStackData to validate
    805                               const pmReadout *output // Output readout
     840                              const pmReadout *output, // Output readout
     841                              const pmReadout *exp    // Exposure map
    806842    )
    807843{
    … …  
    869905    }
    870906
     907    if (exp) {
     908        PM_ASSERT_READOUT_NON_NULL(exp, false);
     909        if (exp->image) {
     910            PS_ASSERT_IMAGES_SIZE_EQUAL(exp->image, output->image, false);
     911        }
     912        if (exp->mask) {
     913            PS_ASSERT_IMAGES_SIZE_EQUAL(exp->mask, output->image, false);
     914        }
     915    }
     916
    871917    return true;
    872918}
    … …  
    937983
    938984/// Constructor
    939 pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance)
     985pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float exp, float addVariance)
    940986{
    941987    pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return
    … …  
    946992    data->inspect = NULL;
    947993    data->weight = weight;
     994    data->exp = exp;
    948995    data->addVariance = addVariance;
    949996
    … …  
    952999
    9531000/// Stack input images
    954 bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
    955                     psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic,
     1001bool pmStackCombine(pmReadout *combined, pmReadout *expmaps, psArray *input,
     1002                    psImageMaskType maskVal, psImageMaskType maskSuspect,
     1003                    psImageMaskType bad, int kernelSize,
     1004                    float iter, float rej, float sys, float olympic,
    9561005                    bool useVariance, bool safe, bool rejection)
    9571006{
    … …  
    9591008    int num;                            // Number of inputs
    9601009    int numCols, numRows;               // Size of (sky) images
    961     if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
     1010    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined, expmaps)) {
    9621011        return false;
    9631012    }
    … …  
    9771026    psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image
    9781027    psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image
     1028    psVector *exps = psVectorAlloc(num, PS_TYPE_F32);    // Exposure times for each image
    9791029    psArray *stack = psArrayAlloc(num); // Stack of readouts
     1030    float totalExpWeight = 0.0;           // Total value of all weighted exposure times
     1031    float totalExp = 0.0;                 // Total exposure time
    9801032    for (int i = 0; i < num; i++) {
    9811033        pmStackData *data = input->data[i]; // Stack data for this input
    9821034        if (!data) {
    9831035            weights->data.F32[i] = 0.0;
     1036            exps->data.F32[i] = NAN;
    9841037            continue;
    9851038        }
    9861039        weights->data.F32[i] = data->weight;
     1040        exps->data.F32[i] = data->exp;
     1041        totalExp += exps->data.F32[i];
     1042        totalExpWeight += exps->data.F32[i] * weights->data.F32[i];
    9871043        pmReadout *ro = data->readout;  // Readout of interest
    9881044        stack->data[i] = psMemIncrRefCounter(ro);
    … …  
    10031059        }
    10041060    }
     1061    totalExpWeight = totalExp / totalExpWeight;    // Convert to inverse
    10051062
    10061063    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
    … …  
    10081065    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
    10091066                                stack)) {
    1010         psError(PS_ERR_UNKNOWN, false, "Input stack is not valid.");
     1067        psError(psErrorCodeLast(), false, "Input stack is not valid.");
    10111068        psFree(stack);
    10121069        return false;
    … …  
    10451102        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    10461103        combinedVariance = combined->variance;
     1104    }
     1105
     1106    psImage *exp = NULL, *expnum = NULL, *expweight = NULL; // Exposure map and exposure number
     1107    if (expmaps) {
     1108        if (!expmaps->image) {
     1109            expmaps->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     1110        }
     1111        exp = expmaps->image;
     1112
     1113        if (!expmaps->mask) {
     1114            expmaps->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     1115        }
     1116        expnum = expmaps->mask;
     1117
     1118        if (!expmaps->variance) {
     1119            expmaps->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     1120        }
     1121        expweight = expmaps->variance;
    10471122    }
    10481123
    … …  
    10831158            bool suspect;               // Suspect pixels in stack?
    10841159            combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
    1085                            input, weights, addVariance, reject, x, y, maskVal, maskSuspect);
    1086             combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe);
     1160                           input, weights, exps, addVariance, reject, x, y, maskVal, maskSuspect);
     1161            combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight,
     1162                          num, buffer, x, y, bad, safe, totalExpWeight);
    10871163
    10881164            if (iter > 0) {
    … …  
    10961172    psFree(weights);
    10971173    psFree(buffer);
    1098 
     1174    psFree(addVariance);
    10991175
    11001176

  • trunk/psModules/src/imcombine/pmStack.h

    r27402 r27427  
    3131    psPixels *inspect;                  ///< Pixels to inspect
    3232    float weight;                       ///< Relative weighting for image
     33    float exp;                          ///< Exposure time
    3334    float addVariance;                  ///< Additional variance when rejecting
    3435} pmStackData;
    … …  
    3738pmStackData *pmStackDataAlloc(pmReadout *readout, ///< Warped readout (sky cell)
    3839                              float weight, ///< Weight to apply
     40                              float exp,    ///< Exposure time
    3941                              float addVariance ///< Additional variance when rejecting
    4042    );
    … …  
    4244/// Stack input images
    4345bool pmStackCombine(pmReadout *combined,///< Combined readout (output)
     46                    pmReadout *expmaps, ///< Exposure maps (output)
    4447                    psArray *input,     ///< Input array of pmStackData
    4548                    psImageMaskType maskVal, ///< Mask value of bad pixels
Note: See TracChangeset for help on using the changeset viewer.