IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34800 for trunk/ppStack/src


Ignore:
Timestamp:
Dec 11, 2012, 2:04:31 PM (14 years ago)
Author:
watersc1
Message:

Merge from my branch including background restoration and stack stage projection cell binned images.

Location:
trunk
Files:
17 edited
3 copied

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/ppStack/src/Makefile.am

    r30620 r34800  
    1 bin_PROGRAMS = ppStack
     1bin_PROGRAMS = ppStack ppStackMedian
    22
    33if HAVE_SVNVERSION
    … …  
    2525ppStack_LDFLAGS = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PSPHOT_LIBS)   $(PPSTATS_LIBS)   $(PPSTACK_LIBS)
    2626
     27ppStackMedian_CFLAGS    = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PPSTATS_CFLAGS) $(PPSTACK_CFLAGS)
     28ppStackMedian_LDFLAGS = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PSPHOT_LIBS)   $(PPSTATS_LIBS)   $(PPSTACK_LIBS)
     29
    2730ppStack_SOURCES =               \
    2831        ppStack.c               \
    … …  
    4346        ppStackCombinePrepare.c \
    4447        ppStackCombineInitial.c \
     48        ppStackCombineAlternate.c       \
     49        ppStackReject.c         \
     50        ppStackCombineFinal.c   \
     51        ppStackCleanup.c        \
     52        ppStackPhotometry.c     \
     53        ppStackFinish.c         \
     54        ppStackTarget.c         \
     55        ppStackUpdateHeader.c   \
     56        ppStackJPEGs.c          \
     57        ppStackStats.c          \
     58        ppStackErrorCodes.c
     59
     60ppStackMedian_SOURCES =         \
     61        ppStackMedian.c         \
     62        ppStackArguments.c      \
     63        ppStackCamera.c         \
     64        ppStackFiles.c          \
     65        ppStackMedianLoop.c             \
     66        ppStackPSF.c            \
     67        ppStackReadout.c        \
     68        ppStackVersion.c        \
     69        ppStackMatch.c          \
     70        ppStackSources.c        \
     71        ppStackThread.c         \
     72        ppStackOptions.c        \
     73        ppStackSetup.c          \
     74        ppStackPrepare.c        \
     75        ppStackConvolve.c       \
     76        ppStackCombinePrepare.c \
     77        ppStackCombineInitial.c \
     78        ppStackCombineAlternate.c       \
    4579        ppStackReject.c         \
    4680        ppStackCombineFinal.c   \

  • trunk/ppStack/src/ppStack.h

    r34092 r34800  
    3838    PPSTACK_FILES_STACK,                // Stack files
    3939    PPSTACK_FILES_UNCONV,               // Unconvolved stack files
    40     PPSTACK_FILES_PHOT                  // Files for photometry
     40    PPSTACK_FILES_PHOT,                 // Files for photometry
     41    PPSTACK_FILES_BKG,                  // Files for bkg
     42    PPSTACK_FILES_MEDIAN_IN,                // Files for median only stacks.
     43    PPSTACK_FILES_MEDIAN_OUT                // Files for median only stacks.
    4144} ppStackFileList;
    4245
    … …  
    111114bool ppStackReadoutFinalThread(psThreadJob *job // Job to process
    112115    );
     116// Perform median stacking for background
     117bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
     118bool ppStackCombineMedian(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
    113119
    114120// Return software version

  • trunk/ppStack/src/ppStackArguments.c

    r31422 r34800  
    324324    }
    325325
     326   
    326327    psTrace("ppStack", 1, "Done parsing arguments\n");
    327328

  • trunk/ppStack/src/ppStackCamera.c

    r30620 r34800  
    156156            psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF
    157157            psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources
    158 
     158            psString bkgmodel = psMetadataLookupStr(&mdok, input, "BKGMODEL"); // Name of warped background model
     159           
    159160            pmFPAfile *imageFile = defineFile(config, NULL, "PPSTACK.INPUT",
    160161                                              image, PM_FPA_FILE_IMAGE); // File for image
    … …  
    215216            }
    216217
     218            // Grab bkgmodel information here
     219            if (!bkgmodel ||  strlen(bkgmodel) == 0) {
     220              // We have no background models.
     221            }
     222            else {
     223              pmFPAfile *inputBKG = defineFile(config,NULL,"PPSTACK.INPUT.BKGMODEL",bkgmodel,
     224                                               PM_FPA_FILE_IMAGE);
     225              if (!inputBKG) {
     226                psError(psErrorCodeLast(), false,
     227                        "Unable to define file from bkgmodel %d (%s)",i,bkgmodel);
     228                return(false);
     229              }
     230            }// End bkgmodel
     231
     232
     233           
     234           
    217235            i++;
    218236        }
    … …  
    460478    jpeg2->save = true;
    461479
     480    // Output background
     481    pmFPAfile *outBkg = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGMODEL");
     482    if (!outBkg) {
     483      psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGMODEL"));
     484      return(false);
     485    }
     486    if (outBkg->type != PM_FPA_FILE_IMAGE) {
     487      psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGMODEL is not of type IMAGE");
     488      return(false);
     489    }
     490    outBkg->save = true;
     491    pmFPAfile *outBkgRest = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGREST");
     492    if (!outBkgRest) {
     493      psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGREST"));
     494      return(false);
     495    }
     496    if (outBkgRest->type != PM_FPA_FILE_IMAGE) {
     497      psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGREST is not of type IMAGE");
     498      return(false);
     499    }
     500   
    462501    // For photometry, we operate on the chip-mosaicked image
    463502    // we create a copy of the mosaicked image for psphot so we can write out a clean image
    … …  
    496535    }
    497536
     537    // Define output file here.
     538   
    498539    return true;
    499540}

  • trunk/ppStack/src/ppStackCleanup.c

    r30874 r34800  
    2222    options->outRO = NULL;
    2323
    24     options->expRO->data_exists = false;
    25     options->expRO->parent->data_exists = false;
    26     options->expRO->parent->parent->data_exists = false;
    27     psFree(options->expRO);
    28     options->expRO = NULL;
     24    if (options->expRO) {
     25      options->expRO->data_exists = false;
     26      if (options->expRO->parent) {
     27        options->expRO->parent->data_exists = false;
     28        options->expRO->parent->parent->data_exists = false;
     29      }
     30      psFree(options->expRO);
     31      options->expRO = NULL;
     32    }
    2933
     34    if (options->bkgRO) {
     35      options->bkgRO->data_exists = false;
     36      if (options->bkgRO->parent) {
     37        options->bkgRO->parent->data_exists = false;
     38        options->bkgRO->parent->parent->data_exists = false;
     39      }
     40      psFree(options->bkgRO);
     41      options->bkgRO = NULL;
     42    }
     43   
    3044    for (int i = 0; i < options->num; i++) {
    3145        pmCellFreeData(options->cells->data[i]);
    … …  
    6276
    6377bool ppStackCleanup (pmConfig *config, ppStackOptions *options) {
    64 
    6578    psExit exitValue = ppStackExitCode(PS_EXIT_SUCCESS); // Exit code
    66 
     79   
    6780    // Ensure everything closes
    6881    if (config) {
    … …  
    7285        ppStackFileActivation(config, PPSTACK_FILES_UNCONV, true);
    7386        ppStackFileActivation(config, PPSTACK_FILES_PHOT, true);
     87        ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_IN, true);
     88        ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_OUT, true);
    7489        if (!ppStackFilesIterateUp(config)) {
    7590            psError(psErrorCodeLast(), false, "Unable to close files.");

  • trunk/ppStack/src/ppStackCombinePrepare.c

    r30620 r34800  
    11#include "ppStack.h"
    22
    3 bool ppStackCombinePrepare(const char *outName, const char *expName,
     3bool ppStackCombinePrepare(const char *outName, const char *expName, const char *bkgName,
    44                           ppStackFileList files, ppStackThreadData *stack,
    55                           ppStackOptions *options, pmConfig *config)
    … …  
    2727    options->outRO = pmReadoutAlloc(cell); // Output readout
    2828
    29     pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell
    30     options->expRO = pmReadoutAlloc(expCell); // Output readout
     29    if (expName) {
     30      pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell
     31      options->expRO = pmReadoutAlloc(expCell); //Output readout
     32    }
     33/*     else { */
     34/*       options->expRO = NULL; */
     35/*     } */
     36    pmCell *bkgCell;
     37    int bkg_r0,bkg_c0;
     38    int bkg_nC,bkg_nR;
     39    if (bkgName) {
     40      ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, true, 0);
     41      pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", 0);
     42      pmFPAview *view = ppStackFilesIterateDown(config);
     43      pmReadout *ro = pmFPAviewThisReadout(view,file->fpa);
     44
     45      bkg_r0 = ro->image->row0;
     46      bkg_c0 = ro->image->col0;
     47      bkg_nC = ro->image->numCols;
     48      bkg_nR = ro->image->numRows;
     49      bkgCell = pmFPAfileThisCell(config->files, view, bkgName); // Bkg cell
     50     
     51      options->bkgRO = pmReadoutAlloc(bkgCell); // BKG readout
     52      //      if (!pmHDUGenerateForFPA(options->bkgRO->parent->parent->parent)) {
     53      options->bkgRO->parent->parent->parent->hdu = pmHDUAlloc(NULL);
     54      if (!options->bkgRO->parent->parent->parent->hdu) {
     55        fprintf(stderr,"failed to generate a HDU for this thing.\n");
     56      }
     57      options->bkgRO->parent->parent->parent->hdu->header = psMetadataCopy(options->bkgRO->parent->parent->parent->hdu->header,
     58                                                                           ro->parent->parent->parent->hdu->header);
     59
     60      options->bkgRO->parent->concepts = psMetadataCopy(options->bkgRO->parent->concepts,
     61                                                        ro->parent->concepts);
     62      options->bkgRO->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->concepts,
     63                                                                ro->parent->parent->concepts);
     64      options->bkgRO->parent->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->parent->concepts,
     65                                                                        ro->parent->parent->parent->concepts);
     66
     67    }
     68    else {
     69      options->bkgRO = NULL;
     70    }
    3171
    3272    psFree(view);
    … …  
    4282    }
    4383
    44     if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) {
     84    if (expName) {
     85      if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) {
    4586        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
    4687        return false;
     88      }
     89    }
     90
     91    if (bkgName) {
     92      if (!pmReadoutStackDefineOutput(options->bkgRO, bkg_c0, bkg_r0, bkg_nC, bkg_nR, false, false, 0)) {
     93        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
     94        return false;
     95      }
     96     
     97
     98     
    4799    }
    48100

  • trunk/ppStack/src/ppStackFiles.c

    r31158 r34800  
    1414/// Files required for the convolution
    1515static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
     16
     17/// Files required for the background
     18static char *filesBkg[] = { "PPSTACK.INPUT.BKGMODEL", NULL };
     19
     20/// Files required for median only stacking
     21static char *filesMedianIn[] =  { "PPSTACK.INPUT",                                                           
     22                              NULL };
     23/// Files required for median only stacking
     24static char *filesMedianOut[] =  { "PPSTACK.OUTPUT",                                                         
     25                              NULL };
    1626
    1727/// Regular (convolved) stack files
    … …  
    1929                              "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT",
    2030                              "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2",
     31                              "PPSTACK.OUTPUT.BKGMODEL",
    2132                              NULL };
    2233/// Unconvolved stack files
    … …  
    4152      case PPSTACK_FILES_UNCONV:   return filesUnconv;
    4253      case PPSTACK_FILES_PHOT:     return filesPhot;
     54    case PPSTACK_FILES_BKG:        return filesBkg;
     55    case PPSTACK_FILES_MEDIAN_IN:        return filesMedianIn;
     56    case PPSTACK_FILES_MEDIAN_OUT:        return filesMedianOut;
    4357      default:
    4458        psAbort("Unrecognised file list: %x", list);

  • trunk/ppStack/src/ppStackLoop.c

    r31435 r34800  
    6363
    6464    // Prepare for combination
    65     if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK, stack, options, config)) {
     65    if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.BKGMODEL", PPSTACK_FILES_STACK, stack, options, config)) {
    6666        psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
    6767        psFree(stack);
    … …  
    8484        pmCellFreeData(options->cells->data[i]);
    8585    }
    86     psFree(stack);
     86    //    psFree(stack);
    8787
    8888    // Pixel rejection
    … …  
    142142    }
    143143
     144    // Generate median background stack here.
     145    if (!ppStackCombineBackground(stack, options, config)) {
     146      psError(psErrorCodeLast(), false, "Unable to generate median of background images.");
     147      psFree(stack);
     148      return false;
     149    }
     150    ppStackFileActivation(config,PPSTACK_FILES_BKG ,false);   
    144151    // Photometry
    145152    psTrace("ppStack", 1, "Photometering stacked image....\n");
    … …  
    171178        return false;
    172179    }
    173     psFree(stack);
     180    //    psFree(stack);
    174181    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
    175182    ppStackMemDump("cleanup");
    … …  
    186193
    187194        // Prepare for combination
    188         if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV,
     195        if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV,
    189196                                   stack, options, config)) {
    190197            psError(psErrorCodeLast(), false, "Unable to prepare for combination.");

  • trunk/ppStack/src/ppStackLoop.h

    r30620 r34800  
    44// Loop over the inputs, doing the combination
    55bool ppStackLoop(
     6    pmConfig *config,                    // Configuration
     7    ppStackOptions *options             // Options for stacking
     8    );
     9// Median only loop.
     10bool ppStackMedianLoop(
    611    pmConfig *config,                    // Configuration
    712    ppStackOptions *options             // Options for stacking
    … …  
    3237    const char *outName,                // Name of output file
    3338    const char *expName,                // Name of exposure file
     39    const char *bkgName,                // Name of background file
    3440    ppStackFileList files,              // Files of interest
    3541    ppStackThreadData *stack,           // Stack

  • trunk/ppStack/src/ppStackOptions.h

    r34093 r34800  
    4848    // Rejection
    4949    psArray *rejected;                  // Rejected pixels
     50    // Background
     51    pmReadout *bkgRO;                   // Output background readout
     52    psArray   *bkgImages;               // Input background images
    5053} ppStackOptions;
    5154

  • trunk/ppStack/src/ppStackPrepare.c

    r34095 r34800  
    132132    psVectorInit(options->inputMask, 0);
    133133    options->exposures = psVectorAlloc(options->num, PS_TYPE_F32);
     134
    134135    psVectorInit(options->exposures, NAN);
    135136
    … …  
    150151
    151152        options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE");
     153        if ((options->exposures->data.F32[i] == 0)||
     154            (!(isfinite(options->exposures->data.F32[i])))){
     155          options->exposures->data.F32[i] = psMetadataLookupF32(NULL,recipe,"DEFAULT.EXPTIME");
     156        }
    152157        options->sumExposure += options->exposures->data.F32[i];
    153158

  • trunk/ppStack/src/ppStackReadout.c

    r30620 r34800  
    194194
    195195
    196 
    197196bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts,
    198197                         const psVector *mask, const psArray *rejected, const psVector *weightings,

  • trunk/ppStack/src/ppStackSetup.c

    r31158 r34800  
    8989    options->origMasks = psArrayAlloc(num);
    9090    options->origVariances = psArrayAlloc(num);
     91    options->bkgImages = psArrayAlloc(num);
    9192    pmFPAview *view = pmFPAviewAlloc(0);
     93    int nullMasks = 0;
     94    int nullVariances = 0;
    9295    for (int i = 0; i < num; i++) {
    9396        {
    … …  
    97100        {
    98101            // We want the convolved mask, since that defines the area that has been tested for outliers
     102          if (options->convolve) {
    99103            options->origMasks->data[i] = psMemIncrRefCounter(options->convMasks->data[i]);
     104          }
     105          else {
     106            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.MASK", i);
     107            options->origMasks->data[i] = pmFPAfileName(file, view, config);
     108          }
     109          if (!(options->origMasks->data[i])) {
     110            nullMasks++;
     111          }
    100112        }
    101113        {
    102114            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i);
    103115            options->origVariances->data[i] = pmFPAfileName(file, view, config);
     116            if (!(options->origVariances->data[i])) {
     117              nullVariances++;
     118            }
    104119        }
     120/*      { */
     121/*        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); */
     122/*        options->bkgImages->data[i] = pmFPAfileName(file, view, config); */
     123/*      } */
    105124    }
     125    if (nullMasks == num) {
     126      psFree(options->origMasks);
     127    }
     128    if (nullVariances == num) {
     129      psFree(options->origVariances);
     130    }
     131   
    106132    psFree(view);
    107133

  • trunk/ppStack/src/ppStackThread.c

    r31158 r34800  
    4040    psFree(stack->maskFits);
    4141    psFree(stack->varianceFits);
     42    psFree(stack->bkgFits);
    4243    return;
    4344}
    … …  
    6768        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
    6869    }
    69 
     70   
    7071    ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
    7172    psMemSetDeallocator(stack, (psFreeFunc)stackThreadDataFree);
    … …  
    7778    stack->maskFits   = psArrayAlloc(numInputs);
    7879    stack->varianceFits = psArrayAlloc(numInputs);
     80    stack->bkgFits    = psArrayAlloc(numInputs);
    7981    for (int i = 0; i < numInputs; i++) {
    8082        if (!cells->data[i]) {
    … …  
    9597            psFree(resolved); \
    9698        }
    97 
     99               
    98100        IMAGE_OPEN(imageNames, stack->imageFits, i);
    99         IMAGE_OPEN(maskNames, stack->maskFits, i);
    100         IMAGE_OPEN(varianceNames, stack->varianceFits, i);
     101        if (maskNames) {
     102          IMAGE_OPEN(maskNames, stack->maskFits, i);
     103        }
     104        if (varianceNames) {
     105          IMAGE_OPEN(varianceNames, stack->varianceFits, i);
     106        }
    101107    }
    102108

  • trunk/ppStack/src/ppStackThread.h

    r31158 r34800  
    2525    psArray *maskFits;                  // FITS file pointers for masks
    2626    psArray *varianceFits;              // FITS file pointers for variances
     27    psArray *bkgFits;                   // FITS file pointers for background models
    2728} ppStackThreadData;
    2829

  • trunk/ppStack/src/ppStackUpdateHeader.c

    r31158 r34800  
    66
    77    pmReadout *outRO = options->outRO;                                      // Output readout
     8    pmReadout *expRO = options->expRO;
    89
    910    // Propagate WCS
    1011    bool wcsDone = false;           // Have we done the WCS?
    1112    for (int i = 0; i < options->num && !wcsDone; i++) {
    12         if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
     13      if ((options->inputMask)&&(options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i])) {
    1314            continue;
    1415        }
    … …  
    123124        snprintf (field, 64, "AIR_%04d", i);
    124125        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image airmass", value);
    125     }   
     126    }
     127
     128    // Copy information into expRO, because it should be there too.
     129    if ((expRO)&&(expRO->parent)) {
     130      pmHDU *expROhdu = pmHDUFromCell(expRO->parent);
     131      if (!expROhdu->header) {
     132        expROhdu->header = psMetadataAlloc();
     133      }
     134      expRO->parent->parent->parent->hdu->header = psMetadataCopy(expRO->parent->parent->parent->hdu->header,
     135                                                                  outRO->parent->parent->parent->hdu->header);
     136     
     137      expRO->parent->concepts = psMetadataCopy(expRO->parent->concepts,
     138                                               outRO->parent->concepts);
     139      expRO->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->concepts,
     140                                                       outRO->parent->concepts);
     141      expRO->parent->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->parent->concepts,
     142                                                               outRO->parent->parent->parent->concepts);
     143    }
     144   
    126145    return true;
    127146}
Note: See TracChangeset for help on using the changeset viewer.