IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27400 for trunk/ppStack/src


Ignore:
Timestamp:
Mar 22, 2010, 8:34:28 PM (16 years ago)
Author:
Paul Price
Message:

Adding exposure map (both exposure time and number of inputs) to the outputs of ppStack. These are calculated both for the regular (convolved) and unconvolved stacks. These should compress pretty well, but I haven't enabled the compression yet. I also plan to add a weighted exposure time output.

Location:
trunk/ppStack/src
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r27319 r27400  
    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/ppStackCamera.c

    r27075 r27400  
    283283    }
    284284
     285
     286    // Exposure image
     287    pmFPA *expFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output
     288    if (!expFPA) {
     289        psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration.");
     290        return false;
     291    }
     292    pmFPAfile *exp = pmFPAfileDefineOutput(config, expFPA, "PPSTACK.OUTPUT.EXP");
     293    psFree(expFPA);                        // Drop reference
     294    if (!exp) {
     295        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.EXP"));
     296        return false;
     297    }
     298    if (exp->type != PM_FPA_FILE_IMAGE) {
     299        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXP is not of type IMAGE");
     300        return false;
     301    }
     302    exp->save = true;
     303
     304    if (!pmFPAAddSourceFromFormat(expFPA, "Stack", exp->format)) {
     305        psError(psErrorCodeLast(), false, "Unable to generate output FPA.");
     306        return false;
     307    }
     308
     309    // Exposure numbers
     310    pmFPAfile *expNum = pmFPAfileDefineOutput(config, exp->fpa, "PPSTACK.OUTPUT.EXPNUM");
     311    if (!expNum) {
     312        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.MASK"));
     313        return false;
     314    }
     315    if (expNum->type != PM_FPA_FILE_MASK) {
     316        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXPNUM is not of type MASK");
     317        return false;
     318    }
     319    expNum->save = true;
     320
     321
     322
    285323    if (havePSFs) {
    286324        pmFPA *psfFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain PSF
    … …  
    302340    }
    303341
    304 #if 1
    305342    // Unconvolved stack
    306343    pmFPA *unconvFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain unconvolved output
    … …  
    351388        unconvVariance->save = true;
    352389    }
    353 #endif
     390
     391
     392    // Exposure image
     393    pmFPA *unconvExpFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output
     394    if (!unconvExpFPA) {
     395        psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration.");
     396        return false;
     397    }
     398    pmFPAfile *unconvExp = pmFPAfileDefineOutput(config, unconvExpFPA, "PPSTACK.UNCONV.EXP");
     399    psFree(unconvExpFPA);               // Drop reference
     400    if (!unconvExp) {
     401        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.EXP"));
     402        return false;
     403    }
     404    if (unconvExp->type != PM_FPA_FILE_IMAGE) {
     405        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXP is not of type IMAGE");
     406        return false;
     407    }
     408    unconvExp->save = true;
     409
     410    if (!pmFPAAddSourceFromFormat(unconvExpFPA, "Stack", unconvExp->format)) {
     411        psError(psErrorCodeLast(), false, "Unable to generate output FPA.");
     412        return false;
     413    }
     414
     415    // Output mask
     416    pmFPAfile *unconvExpNum = pmFPAfileDefineOutput(config, unconvExp->fpa, "PPSTACK.UNCONV.EXPNUM");
     417    if (!unconvExpNum) {
     418        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.MASK"));
     419        return false;
     420    }
     421    if (unconvExpNum->type != PM_FPA_FILE_MASK) {
     422        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXPNUM is not of type MASK");
     423        return false;
     424    }
     425    unconvExpNum->save = true;
     426
    354427
    355428    // Output JPEGs

  • trunk/ppStack/src/ppStackCleanup.c

    r27343 r27400  
    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

    r27319 r27400  
    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
    … …  
    4345    }
    4446
    45     if (!target->mask) {
    46         target->mask = psImageAlloc(target->image->numCols, target->image->numRows, PS_TYPE_IMAGE_MASK);
     47    if (!outRO->mask) {
     48        outRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     49    }
     50    if (!expRO->mask) {
     51        expRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    4752    }
    4853
    … …  
    6267        }
    6368
    64         // call: ppStackReadoutFinal(config, target, readouts, rejected)
     69        // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
    6570        psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
    66         psArrayAdd(job->args, 1, target);
    6771        psArrayAdd(job->args, 1, thread);
    6872        psArrayAdd(job->args, 1, reject);
    … …  
    105109        }
    106110        if (sumWeights > 0.0) {
    107             target->covariance = psImageCovarianceSum(covariances);
    108             psBinaryOp(target->covariance->image, target->covariance->image, "/",
     111            outRO->covariance = psImageCovarianceSum(covariances);
     112            psBinaryOp(outRO->covariance->image, outRO->covariance->image, "/",
    109113                       psScalarAlloc(sumWeights, PS_TYPE_F32));
    110             psImageCovarianceTransfer(target->variance, target->covariance);
     114            psImageCovarianceTransfer(outRO->variance, outRO->covariance);
    111115        }
    112116    } else {
    113         target->covariance = psImageCovarianceNone();
     117        outRO->covariance = psImageCovarianceNone();
    114118    }
    115119
    … …  
    127131            wcsDone = true;
    128132            pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
    129             pmHDU *outHDU = pmHDUFromCell(target->parent); // Output HDU
    130             pmChip *outChip = target->parent->parent; // Output chip
     133            pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
     134            pmChip *outChip = outRO->parent->parent; // Output chip
    131135            pmFPA *outFPA = outChip->parent; // Output FPA
    132136            if (!outHDU || !inHDU) {
    … …  
    151155    }
    152156
     157    // Set exposure time correctly
     158    {
     159        float exptime = 0.0;            // Summed exposure time
     160        for (int i = 0; i < options->num; i++) {
     161            if (options->inputMask) {
     162                continue;
     163            }
     164            exptime += options->exposures->data.F32[i];
     165        }
     166
     167        {
     168            psMetadataItem *item = psMetadataLookup(outRO->parent->concepts, "CELL.EXPOSURE");
     169            item->data.F32 = exptime;
     170        }
     171        {
     172            psMetadataItem *item = psMetadataLookup(outRO->parent->parent->parent->concepts, "FPA.EXPOSURE");
     173            item->data.F32 = exptime;
     174        }
     175    }
     176
    153177    // Put version information into the header
    154     pmHDU *hdu = pmHDUFromCell(target->parent);
     178    pmHDU *hdu = pmHDUFromCell(outRO->parent);
    155179    if (!hdu) {
    156180        psError(PPSTACK_ERR_PROG, false, "Unable to find HDU for output.");
    … …  
    168192    psStringAppend(&name, "combined_image_final_%d.fits", pass);
    169193    pass++;
    170     ppStackWriteImage(name, NULL, target->image, config);
     194    ppStackWriteImage(name, NULL, outRO->image, config);
    171195    psStringSubstitute(&name, "mask", "image");
    172     ppStackWriteImage(name, NULL, target->mask, config);
     196    ppStackWriteImage(name, NULL, outRO->mask, config);
    173197    psStringSubstitute(&name, "variance", "mask");
    174     ppStackWriteImage(name, NULL, target->variance, config);
     198    ppStackWriteImage(name, NULL, outRO->variance, config);
    175199    psFree(name);
    176200
    177     pmStackVisualPlotTestImage(target->image, "combined_image_final.fits");
     201    pmStackVisualPlotTestImage(outRO->image, "combined_image_final.fits");
    178202#endif
    179203

  • trunk/ppStack/src/ppStackCombineInitial.c

    r27309 r27400  
    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

    r27319 r27400  
    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/ppStackFiles.c

    r27319 r27400  
    2323/// Regular (convolved) stack files
    2424static char *filesStack[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
    25                               "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
     25                              "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM",
     26                              "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2",
     27                              NULL };
    2628/// Unconvolved stack files
    27 static char *filesUnconv[] = { "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE", NULL };
     29static char *filesUnconv[] = { "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE",
     30                               "PPSTACK.UNCONV.EXP", "PPSTACK.UNCONV.EXPNUM", NULL };
    2831
    2932/// Files for photometry

  • trunk/ppStack/src/ppStackLoop.c

    r27343 r27400  
    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

    r27319 r27400  
    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

    r27218 r27400  
    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

    r27218 r27400  
    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

    r27075 r27400  
    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

    r27319 r27400  
    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

    r27319 r27400  
    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

    r27329 r27400  
    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

    r27319 r27400  
    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);
Note: See TracChangeset for help on using the changeset viewer.