IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 17, 2009, 12:08:50 PM (17 years ago)
Author:
beaumont
Message:

merged with head

Location:
branches/cnb_branches/cnb_branch_20090301
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branches/cnb_branch_20090301

  • branches/cnb_branches/cnb_branch_20090301/ppStack

  • branches/cnb_branches/cnb_branch_20090301/ppStack/src/ppStackLoop.c

    r21526 r23352  
    44
    55#include <stdio.h>
    6 #include <unistd.h>
    7 #include <string.h>
    8 #include <libgen.h>
    96#include <pslib.h>
    107#include <psmodules.h>
    11 #include <ppStats.h>
    128
    139#include "ppStack.h"
    14 
    15 //#define TESTING
    16 
    17 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
    18 #define PIXELS_BUFFER 1024              // Initial size of pixel lists
    19 
    20 // Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
    21 
    22 // Files required in preparation for convolution
    23 static char *prepareFiles[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL };
    24 
    25 // Files required for the convolution
    26 static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
    27 
    28 // Output files for the combination
    29 static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
    30                                 "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
    31 
    32 // Files for photometry
    33 static char *photFiles[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
    34                              "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
    35                              "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
    36                              "PSPHOT.INPUT.CMF", NULL };
    37 
    38 static void memDump(const char *name)
    39 {
    40     return;
    41 
    42     static int num = 0;                 // Counter, to make files unique and give an idea of sequence
    43 
    44     psString filename = NULL;           // Name of file
    45     psStringAppend(&filename, "memdump_%s_%03d.txt", name, num);
    46     FILE *memFile = fopen(filename, "w");
    47     psFree(filename);
    48 
    49     psMemBlock **leaks = NULL;
    50     int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
    51     fprintf(memFile, "# MemBlock Size Source\n");
    52     unsigned long total = 0;            // Total memory used
    53     for (int i = 0; i < numLeaks; i++) {
    54         psMemBlock *mb = leaks[i];
    55         fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
    56                 mb->file, mb->lineno);
    57         total += mb->userMemorySize;
    58     }
    59     fclose(memFile);
    60     psFree(leaks);
    61 
    62     fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, sbrk(0));
    63     num++;
    64 }
    65 
    66 
    67 
    68 // Activate/deactivate a list of files
    69 static void fileActivation(pmConfig *config, // Configuration
    70                            char **files, // Files to turn on/off
    71                            bool state   // Activation state
    72     )
    73 {
    74     assert(config);
    75     assert(files);
    76 
    77     for (int i = 0; files[i] != NULL; i++) {
    78         pmFPAfileActivate(config->files, state, files[i]);
    79     }
    80     return;
    81 }
    82 
    83 // Activate/deactivate a single element for a list
    84 static void fileActivationSingle(pmConfig *config, // Configuration
    85                                  char **files, // Files to turn on/off
    86                                  bool state,   // Activation state
    87                                  int num // Number of file in sequence
    88                                  )
    89 {
    90     assert(config);
    91     assert(files);
    92 
    93     for (int i = 0; files[i] != NULL; i++) {
    94         pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
    95     }
    96     return;
    97 }
    98 
    99 // Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
    100 static pmFPAview *filesIterateDown(pmConfig *config // Configuration
    101                                   )
    102 {
    103     assert(config);
    104 
    105     pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
    106     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    107         return NULL;
    108     }
    109     view->chip = 0;
    110     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    111         return NULL;
    112     }
    113     view->cell = 0;
    114     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    115         return NULL;
    116     }
    117     view->readout = 0;
    118     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    119         return NULL;
    120     }
    121     return view;
    122 }
    123 
    124 // Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells
    125 static bool filesIterateUp(pmConfig *config // Configuration
    126                            )
    127 {
    128     assert(config);
    129 
    130     pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
    131     view->chip = view->cell = view->readout = 0;
    132     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    133         return false;
    134     }
    135     view->readout = -1;
    136     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    137         return false;
    138     }
    139     view->cell = -1;
    140     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    141         return false;
    142     }
    143     view->chip = -1;
    144     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    145         return false;
    146     }
    147     psFree(view);
    148     return true;
    149 }
    150 
    151 // Write an image to a FITS file
    152 static bool writeImage(const char *name, // Name of image
    153                        psMetadata *header, // Header
    154                        const psImage *image, // Image
    155                        pmConfig *config // Configuration
    156                        )
    157 {
    158     assert(name);
    159     assert(image);
    160 
    161     psString resolved = pmConfigConvertFilename(name, config, true, true); // Resolved file name
    162     psFits *fits = psFitsOpen(resolved, "w");
    163     if (!fits) {
    164         psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved);
    165         psFree(resolved);
    166         return false;
    167     }
    168     if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
    169         psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved);
    170         psFitsClose(fits);
    171         psFree(resolved);
    172         return false;
    173     }
    174     psFitsClose(fits);
    175     psFree(resolved);
    176     return true;
    177 }
    178 
     10#include "ppStackLoop.h"
    17911
    18012bool ppStackLoop(pmConfig *config)
    … …  
    18315
    18416    psTimerStart("PPSTACK_TOTAL");
     17    psTimerStart("PPSTACK_STEPS");
    18518
    186     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    187     psAssert(recipe, "We've thrown an error on this before.");
     19    ppStackOptions *options = ppStackOptionsAlloc(); // Options for stacking
    18820
    189     bool mdok;                          // Status of MD lookup
    190     bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files?
    191     const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
    192     if (!tempDir) {
    193         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
     21    // Setup
     22    psTrace("ppStack", 1, "Setup....\n");
     23    if (!ppStackSetup(options, config)) {
     24        psError(PS_ERR_UNKNOWN, false, "Unable to setup.");
     25        psFree(options);
     26        return false;
     27    }
     28    psLogMsg("ppStack", PS_LOG_INFO, "Stage 0: Setup: %f sec", psTimerClear("PPSTACK_STEPS"));
     29    ppStackMemDump("setup");
     30
     31
     32    // Preparation for stacking
     33    psTrace("ppStack", 1, "Preparation for stacking: merging sources, determining target PSF....\n");
     34    if (!ppStackPrepare(options, config)) {
     35        psError(PS_ERR_UNKNOWN, false, "Unable to prepare for stacking.");
     36        psFree(options);
     37        return false;
     38    }
     39    psLogMsg("ppStack", PS_LOG_INFO, "Stage 1: Load Sources and Generate Target PSF: %f sec",
     40             psTimerClear("PPSTACK_STEPS"));
     41    ppStackMemDump("prepare");
     42
     43
     44    // Convolve inputs
     45    psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
     46    if (!ppStackConvolve(options, config)) {
     47        psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
     48        psFree(options);
     49        return false;
     50    }
     51    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec",
     52             psTimerClear("PPSTACK_STEPS"));
     53    ppStackMemDump("convolve");
     54
     55
     56    // Start threading
     57    ppStackThreadInit();
     58    ppStackThreadData *stack = ppStackThreadDataSetup(options, config);
     59    if (!stack) {
     60        psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
     61        psFree(options);
    19462        return false;
    19563    }
    19664
    197     psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, "OUTPUT")); // Name for temporary files
    198     const char *tempName = basename(outputName);
    199     if (!tempName) {
    200         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
     65    // Initial combination
     66    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
     67    if (!ppStackCombineInitial(stack, options, config)) {
     68        psError(PS_ERR_UNKNOWN, false, "Unable to perform initial combination.");
     69        psFree(stack);
     70        psFree(options);
    20171        return false;
    20272    }
     73    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec",
     74             psTimerClear("PPSTACK_STEPS"));
     75    ppStackMemDump("convolve");
     76    psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
    20377
    204     const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images
    205     const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks
    206     const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp variance maps
    207     if (!tempImage || !tempMask || !tempVariance) {
    208         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    209                 "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
     78
     79    // Pixel rejection
     80    psTrace("ppStack", 1, "Reject pixels....\n");
     81    if (!ppStackReject(options, config)) {
     82        psError(PS_ERR_UNKNOWN, false, "Unable to reject pixels.");
     83        psFree(stack);
     84        psFree(options);
    21085        return false;
    21186    }
     87    psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS"));
     88    ppStackMemDump("reject");
    21289
    213     float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution
    214     float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject
    215                                                                      // before rejecting entire image
    216     float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
    217 
    218     const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics
    219     psMetadata *stats = NULL;           // Container for statistics
    220     FILE *statsFile = NULL;             // File stream for statistics
    221     if (statsName && strlen(statsName) > 0) {
    222         psString resolved = pmConfigConvertFilename(statsName, config, true, true); // Resolved filename
    223         statsFile = fopen(resolved, "w");
    224         if (!statsFile) {
    225             psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
    226             psFree(resolved);
    227             return false;
    228         } else {
    229             stats = psMetadataAlloc();
    230         }
    231         psFree(resolved);
    232     }
    233 
    234     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
    235     if (!output) {
    236         psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!");
    237         return false;
    238     }
    239     int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
    240 
    241     psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    242     int overlap = 2 * psMetadataLookupS32(NULL, ppsub,
    243                                           "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
    244 
    245     if (!pmConfigMaskSetBits(NULL, NULL, config)) {
    246         psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
    247         return false;
    248     }
    249 
    250     memDump("start");
    251 
    252     psLogMsg("ppStack", PS_LOG_INFO, "Stage 0 : Initialization and Configuration : %f sec", psTimerClear("PPSTACK_STEPS"));
    253 
    254     // Preparation iteration: Load the sources, and get a target PSF model
    255     psTrace("ppStack", 1, "Determining target PSF....\n");
    256     psArray *sourceLists = psArrayAlloc(num); // Individual lists of sources for matching
    257     pmPSF *targetPSF = NULL;            // Target PSF
    258     float sumExposure = NAN;            // Sum of exposure times
    259     if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
    260         pmFPAfileActivate(config->files, false, NULL);
    261         fileActivation(config, prepareFiles, true);
    262         pmFPAview *view = filesIterateDown(config);
    263         if (!view) {
    264             return false;
    265         }
    266 
    267         // Generate list of PSFs
    268         psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
    269                                                                "^PPSTACK.INPUT$"); // Iterator
    270         psMetadataItem *fileItem; // Item from iteration
    271         psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
    272         int numCols = 0, numRows = 0;   // Size of image
    273         int index = 0;                  // Index for file
    274         while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
    275             assert(fileItem->type == PS_DATA_UNKNOWN);
    276             pmFPAfile *inputFile = fileItem->data.V; // An input file
    277             pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
    278             pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
    279             if (!psf) {
    280                 psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
    281                 psFree(view);
    282                 psFree(sourceLists);
    283                 psFree(fileIter);
    284                 psFree(psfs);
    285                 return false;
    286             }
    287             psfs->data[index] = psMemIncrRefCounter(psf);
    288             psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF");
    289 
    290             pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
    291             pmHDU *hdu = pmHDUFromCell(cell);
    292             assert(hdu && hdu->header);
    293             int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns
    294             int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
    295             if (naxis1 <= 0 || naxis2 <= 0) {
    296                 psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
    297                 psFree(view);
    298                 psFree(sourceLists);
    299                 psFree(fileIter);
    300                 psFree(psfs);
    301                 return false;
    302             }
    303             if (numCols == 0 && numRows == 0) {
    304                 numCols = naxis1;
    305                 numRows = naxis2;
    306             }
    307 
    308             pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
    309             psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
    310             if (!sources) {
    311                 psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
    312                 return NULL;
    313             }
    314             sourceLists->data[index] = psMemIncrRefCounter(sources);
    315 
    316             // Re-do photometry if we don't trust the source lists
    317             if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
    318                 psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
    319                 pmFPAfileActivate(config->files, false, NULL);
    320                 fileActivationSingle(config, convolveFiles, true, index);
    321                 pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
    322                 pmFPAview *view = filesIterateDown(config);
    323                 if (!view) {
    324                     psFree(sourceLists);
    325                     psFree(targetPSF);
    326                     return false;
    327                 }
    328 
    329                 pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
    330 
    331                 if (!ppStackInputPhotometry(ro, sources, config)) {
    332                     psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
    333                     psFree(sourceLists);
    334                     psFree(targetPSF);
    335                     return false;
    336                 }
    337 
    338                 psFree(view);
    339                 if (!filesIterateUp(config)) {
    340                     psFree(sourceLists);
    341                     psFree(targetPSF);
    342                     return false;
    343                 }
    344                 pmFPAfileActivate(config->files, false, NULL);
    345                 fileActivation(config, prepareFiles, true);
    346             }
    347 
    348             index++;
    349         }
    350         psFree(fileIter);
    351 
    352         // Zero point calibration
    353         sumExposure = ppStackSourcesTransparency(sourceLists, view, config);
    354         if (!isfinite(sumExposure) || sumExposure <= 0) {
    355             psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
    356             psFree(sourceLists);
    357             psFree(targetPSF);
    358             return false;
    359         }
    360 
    361         // Generate target PSF
    362         targetPSF = ppStackPSF(config, numCols, numRows, psfs);
    363         psFree(psfs);
    364         if (!targetPSF) {
    365             psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
    366             psFree(sourceLists);
    367             psFree(view);
    368             return false;
    369         }
    370         psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
    371                          "Target PSF for stack", targetPSF);
    372 
    373         pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Output chip
    374         psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
    375                          "Target PSF", targetPSF);
    376         outChip->data_exists = true;
    377 
    378         psFree(view);
    379         if (!filesIterateUp(config)) {
    380             return false;
    381         }
    382     }
    383 
    384     psLogMsg("ppStack", PS_LOG_INFO, "Stage 1 : Load Sources and Generate Target PSF : %f sec", psTimerClear("PPSTACK_STEPS"));
    385     memDump("psf");
    386 
    387     psArray *imageNames = psArrayAlloc(num);
    388     psArray *maskNames = psArrayAlloc(num);
    389     psArray *varianceNames = psArrayAlloc(num);
    390     for (int i = 0; i < num; i++) {
    391         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
    392         psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
    393         psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
    394         psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
    395         psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
    396         imageNames->data[i] = imageName;
    397         maskNames->data[i] = maskName;
    398         varianceNames->data[i] = varianceName;
    399     }
    400     // Free the outputName that we grabbed above to set up tempName.
    401     psFree(outputName);
    402 
    403     // Generate convolutions and write them to disk
    404     psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
    405     psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
    406     psArray *subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
    407     psArray *subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking
    408     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    409     int numGood = 0;                    // Number of good frames
    410     int numCols = 0, numRows = 0;       // Size of image
    411     psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
    412     psVectorInit(inputMask, 0);
    413     psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
    414     psVectorInit(matchChi2, NAN);
    415     psVector *weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
    416     psVectorInit(weightings, NAN);
    417     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
    418     psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
    419     psArray *covariances = psArrayAlloc(num); // Covariance matrices
    420     for (int i = 0; i < num; i++) {
    421         psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
    422         pmFPAfileActivate(config->files, false, NULL);
    423         fileActivationSingle(config, convolveFiles, true, i);
    424         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
    425         pmFPAview *view = filesIterateDown(config);
    426         if (!view) {
    427             psFree(sourceLists);
    428             psFree(targetPSF);
    429             psFree(rng);
    430             psFree(inputMask);
    431             psFree(matchChi2);
    432             psFree(fpaList);
    433             psFree(cellList);
    434             psFree(covariances);
    435             return false;
    436         }
    437 
    438         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout
    439         psFree(view);
    440 
    441         if (numCols == 0 && numRows == 0) {
    442             numCols = readout->image->numCols;
    443             numRows = readout->image->numRows;
    444         } else if (numCols != readout->image->numCols || numRows != readout->image->numRows) {
    445             psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d",
    446                     readout->image->numCols, readout->image->numRows, numCols, numRows);
    447             psFree(sourceLists);
    448             psFree(targetPSF);
    449             psFree(rng);
    450             psFree(inputMask);
    451             psFree(matchChi2);
    452             psFree(fpaList);
    453             psFree(cellList);
    454             psFree(covariances);
    455             return false;
    456         }
    457 
    458         // Background subtraction, scaling and normalisation is performed automatically by the image matching
    459         psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
    460         psTimerStart("PPSTACK_MATCH");
    461 
    462         if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
    463                           sourceLists->data[i], targetPSF, rng, config)) {
    464             psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
    465             inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH;
    466             psErrorClear();
    467             continue;
    468         }
    469         covariances->data[i] = psMemIncrRefCounter(readout->covariance);
    470 
    471         if (stats) {
    472             pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
    473                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel
    474             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK,
    475                              "Time to match PSF", psTimerMark("PPSTACK_MATCH"));
    476             psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
    477                              "Mean deviation for stamps", kernels->mean);
    478             psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
    479                              "RMS deviation for stamps", kernels->rms);
    480             psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
    481                              "Number of stamps", kernels->numStamps);
    482             float deconv = psMetadataLookupF32(NULL, readout->analysis,
    483                                                PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction
    484             psMetadataAddF32(stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK,
    485                              "Deconvolution fraction for kernel", deconv);
    486             psMetadataAddF32(stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK,
    487                              "Weighting for image", weightings->data.F32[i]);
    488         }
    489         psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH"));
    490 
    491         subRegions->data[i] = regions;
    492         subKernels->data[i] = kernels;
    493 
    494         // Write the temporary convolved files
    495         pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
    496         assert(hdu);
    497         writeImage(imageNames->data[i], hdu->header, readout->image, config);
    498         psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
    499         pmConfigMaskWriteHeader(config, maskHeader);
    500         writeImage(maskNames->data[i], maskHeader, readout->mask, config);
    501         psFree(maskHeader);
    502         psImageCovarianceTransfer(readout->variance, readout->covariance);
    503         writeImage(varianceNames->data[i], hdu->header, readout->variance, config);
    504 #ifdef TESTING
    505         {
    506             psString name = NULL;
    507             psStringAppend(&name, "covariance_%d.fits", i);
    508             writeImage(name, hdu->header, readout->covariance->image, config);
    509             psFree(name);
    510         }
    511 #endif
    512 
    513         pmCell *inCell = readout->parent; // Input cell
    514 
    515         psListAdd(cellList, PS_LIST_TAIL, inCell);
    516         psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent);
    517 
    518         cells->data[i] = psMemIncrRefCounter(inCell);
    519         if (!filesIterateUp(config)) {
    520             psFree(fpaList);
    521             psFree(cellList);
    522             psFree(covariances);
    523             return false;
    524         }
    525         numGood++;
    526 
    527         memDump("match");
    528     }
    529     psFree(sourceLists);
    530     psFree(targetPSF);
    531     psFree(rng);
    532 
    533     psLogMsg("ppStack", PS_LOG_INFO, "Stage 2 : Generate Convolutions and Save : %f sec", psTimerClear("PPSTACK_STEPS"));
    534 
    535     if (numGood == 0) {
    536         psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
    537         psFree(subKernels);
    538         psFree(subRegions);
    539         psFree(cells);
    540         psFree(inputMask);
    541         psFree(matchChi2);
    542         psFree(fpaList);
    543         psFree(cellList);
    544         psFree(covariances);
    545         return false;
    546     }
    547 
    548 
    549     // Update concepts for output
    550     {
    551         pmFPAview view;                 // View for output
    552         view.chip = view.cell = view.readout = 0;
    553         pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell
    554         pmFPA *outFPA = outCell->parent->parent; // Output FPA
    555         pmConceptsAverageFPAs(outFPA, fpaList);
    556         pmConceptsAverageCells(outCell, cellList, NULL, NULL, true);
    557         psFree(fpaList);
    558         psFree(cellList);
    559 
    560         // Update the value of a concept
    561 #define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \
    562             psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
    563             psAssert(item, "Concept should be present"); \
    564             psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \
    565             item->data.F32 = VALUE; \
    566         }
    567 
    568         UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  sumExposure);
    569         UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", sumExposure);
    570         UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
    571     }
    572 
    573 
    574     // Reject images out-of-hand on the basis of their match chi^2
    575     {
    576         psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort
    577         for (int i = 0; i < num; i++) {
    578             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
    579                 continue;
    580             }
    581             values->data.F32[values->n++] = matchChi2->data.F32[i];
    582         }
    583         assert(values->n == numGood);
    584         if (!psVectorSortInPlace(values)) {
    585             psError(PS_ERR_UNKNOWN, false, "Unable to sort vector.");
    586             psFree(subKernels);
    587             psFree(subRegions);
    588             psFree(cells);
    589             psFree(inputMask);
    590             psFree(matchChi2);
    591             psFree(values);
    592             psFree(covariances);
    593             return false;
    594         }
    595         float median = numGood % 2 ? values->data.F32[numGood / 2] :
    596             0.5 * (values->data.F32[numGood / 2 - 1] + values->data.F32[numGood / 2]);
    597 
    598         float rms = 0.74 * (values->data.F32[numGood * 3 / 4] -
    599                             values->data.F32[numGood / 4]); // Estimated RMS from interquartile range
    600 
    601         psFree(values);
    602 
    603         float rej = psMetadataLookupF32(NULL, recipe, "MATCH.REJ"); // Rejection threshold (stdevs)
    604         if (isfinite(rej)) {
    605             float thresh = median + rej * rms; // Threshold for rejection
    606             psLogMsg("ppStack", PS_LOG_INFO, "chi^2 rejection threshold = %f + %f * %f = %f",
    607                      median, rej, rms, thresh);
    608 
    609             int numRej = 0;                 // Number rejected
    610             numGood = 0;                    // Number of good images
    611             for (int i = 0; i < num; i++) {
    612               if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
    613                     continue;
    614                 }
    615                 if (matchChi2->data.F32[i] > thresh) {
    616                     numRej++;
    617                     inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2;
    618                     psLogMsg("ppStack", PS_LOG_INFO, "Rejecting image %d because of large matching chi^2: %f",
    619                              i, matchChi2->data.F32[i]);
    620                 } else {
    621                     psLogMsg("ppStack", PS_LOG_INFO, "Image %d has matching chi^2: %f",
    622                              i, matchChi2->data.F32[i]);
    623                     numGood++;
    624                 }
    625             }
    626         }
    627     }
    628 
    629     if (numGood == 0) {
    630         psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
    631         psFree(subKernels);
    632         psFree(subRegions);
    633         psFree(cells);
    634         psFree(inputMask);
    635         psFree(matchChi2);
    636         psFree(covariances);
    637         return false;
    638     }
    639 
    640 #ifdef TESTING
    641     psTraceSetLevel("psModules.imcombine", 7);
    642 #endif
    643 
    644     psLogMsg("ppStack", PS_LOG_INFO, "Stage 3 : Basic Image Rejection  : %f sec", psTimerClear("PPSTACK_STEPS"));
    645 
    646     // Start threading
    647     ppStackThreadInit();
    648     ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames,
    649                                                       covariances, config);
    650     if (!stack) {
    651         psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
    652         psFree(subKernels);
    653         psFree(subRegions);
    654         psFree(inputMask);
    655         psFree(matchChi2);
    656         psFree(cells);
    657         psFree(covariances);
    658         return false;
    659     }
    660 
    661     psTimerStart("PPSTACK_INITIAL");
    662 
    663     memDump("preinitial");
    664 
    665     // Stack the convolved files
    666     psTrace("ppStack", 1, "Initial stack of convolved images....\n");
    667     pmReadout *outRO = NULL;            // Output readout
    668     pmFPAview *view = NULL;             // View to readout
    669     psArray *inspect = NULL;            // Array of arrays of pixels to inspect
    670     {
    671         int row0, col0;                 // Offset for readout
    672         int numCols, numRows;           // Size of readout
    673         ppStackThread *thread = stack->threads->data[0]; // Representative thread
    674         if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
    675             psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
    676             psFree(subKernels);
    677             psFree(subRegions);
    678             psFree(stack);
    679             psFree(inputMask);
    680             psFree(matchChi2);
    681             psFree(cells);
    682             psFree(covariances);
    683             return false;
    684         }
    685 
    686         pmFPAfileActivate(config->files, false, NULL);
    687         fileActivation(config, combineFiles, true);
    688         view = filesIterateDown(config);
    689         if (!view) {
    690             psFree(subKernels);
    691             psFree(subRegions);
    692             psFree(stack);
    693             psFree(inputMask);
    694             psFree(matchChi2);
    695             psFree(cells);
    696             psFree(covariances);
    697             return false;
    698         }
    699 
    700         pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
    701         outRO = pmReadoutAlloc(outCell); // Output readout
    702 
    703         psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    704         psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    705         if (!pmReadoutStackDefineOutput(outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
    706             psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
    707             psFree(subKernels);
    708             psFree(subRegions);
    709             psFree(stack);
    710             psFree(inputMask);
    711             psFree(matchChi2);
    712             psFree(view);
    713             psFree(outRO);
    714             psFree(cells);
    715             psFree(covariances);
    716             return false;
    717         }
    718 
    719         psFree(cells);
    720 
    721         bool status;                    // Status of read
    722         int numChunk;                   // Number of chunks
    723         for (numChunk = 0; true; numChunk++) {
    724             ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap);
    725             if (!status) {
    726                 // Something went wrong
    727                 psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
    728                 psFree(subKernels);
    729                 psFree(subRegions);
    730                 psFree(stack);
    731                 psFree(inputMask);
    732                 psFree(matchChi2);
    733                 psFree(view);
    734                 psFree(outRO);
    735                 psFree(covariances);
    736                 return false;
    737             }
    738             if (!thread) {
    739                 // Nothing more to read
    740                 break;
    741             }
    742 
    743             // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
    744             psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
    745             psArrayAdd(job->args, 1, thread);
    746             psArrayAdd(job->args, 1, config);
    747             psArrayAdd(job->args, 1, outRO);
    748             psArrayAdd(job->args, 1, inputMask);
    749             psArrayAdd(job->args, 1, weightings);
    750             psArrayAdd(job->args, 1, matchChi2);
    751             if (!psThreadJobAddPending(job)) {
    752                 psFree(job);
    753                 psFree(subKernels);
    754                 psFree(subRegions);
    755                 psFree(stack);
    756                 psFree(inputMask);
    757                 psFree(weightings);
    758                 psFree(matchChi2);
    759                 psFree(view);
    760                 psFree(outRO);
    761                 psFree(covariances);
    762                 return false;
    763             }
    764             psFree(job);
    765         }
    766 
    767         if (!psThreadPoolWait(false)) {
    768             psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination.");
    769             psFree(subKernels);
    770             psFree(subRegions);
    771             psFree(stack);
    772             psFree(inputMask);
    773             psFree(matchChi2);
    774             psFree(view);
    775             psFree(outRO);
    776             psFree(covariances);
    777             return false;
    778         }
    779 
    780         // Harvest the jobs, gathering the inspection lists
    781         inspect = psArrayAlloc(num);
    782         for (int i = 0; i < num; i++) {
    783             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    784                 continue;
    785             }
    786             inspect->data[i] = psArrayAllocEmpty(numChunk);
    787         }
    788         psThreadJob *job;               // Completed job
    789         while ((job = psThreadJobGetDone())) {
    790             psAssert(strcmp(job->type, "PPSTACK_INITIAL_COMBINE") == 0,
    791                      "Job has incorrect type: %s", job->type);
    792             psArray *results = job->results; // Results of job
    793             for (int i = 0; i < num; i++) {
    794                 if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    795                     continue;
    796                 }
    797                 inspect->data[i] = psArrayAdd(inspect->data[i], 1, results->data[i]);
    798             }
    799             psFree(job);
    800         }
    801 
    802         memDump("initial");
    803     }
    804 
    805 #ifdef TESTING
    806     writeImage("combined_initial.fits", NULL, outRO->image, config);
    807 #endif
    808 
    809     if (stats) {
    810         psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_INITIAL", 0,
    811                          "Time to make initial stack", psTimerMark("PPSTACK_INITIAL"));
    812     }
    813     psLogMsg("ppStack", PS_LOG_INFO, "Stage 4 : Make Initial Stack : %f sec", psTimerClear("PPSTACK_STEPS"));
    814     psLogMsg("ppStack", PS_LOG_INFO, "Time to make initial stack: %f sec", psTimerClear("PPSTACK_INITIAL"));
    815 
    816     psTimerStart("PPSTACK_REJECT");
    817 
    818     // Pixel rejection
    819     psArray *rejected = psArrayAlloc(num);
    820     {
    821         int numRejected = 0;        // Number of inputs rejected completely
    822 
    823         // Count images rejected out of hand
    824         for (int i = 0; i < num; i++) {
    825             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    826                 numRejected++;
    827             }
    828         }
    829 
    830         for (int i = 0; i < num; i++) {
    831             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    832                 continue;
    833             }
    834 
    835             psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start
    836             psArrayAdd(job->args, 1, inspect);
    837             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
    838             if (!psThreadJobAddPending(job)) {
    839                 psFree(job);
    840                 psFree(subKernels);
    841                 psFree(subRegions);
    842                 psFree(stack);
    843                 psFree(inputMask);
    844                 psFree(view);
    845                 psFree(outRO);
    846                 psFree(inspect);
    847                 psFree(rejected);
    848                 psFree(covariances);
    849                 psFree(matchChi2);
    850                 return false;
    851             }
    852             psFree(job);
    853         }
    854 
    855         if (!psThreadPoolWait(true)) {
    856             psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists.");
    857             psFree(subKernels);
    858             psFree(subRegions);
    859             psFree(stack);
    860             psFree(inputMask);
    861             psFree(view);
    862             psFree(outRO);
    863             psFree(inspect);
    864             psFree(rejected);
    865             psFree(covariances);
    866             psFree(matchChi2);
    867             return false;
    868         }
    869 
    870         if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) {
    871             pmStackRejectThreadsInit();
    872         }
    873 
    874         int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches
    875 
    876         // Reject bad pixels
    877         for (int i = 0; i < num; i++) {
    878             if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    879                 continue;
    880             }
    881             psTimerStart("PPSTACK_REJECT");
    882 
    883 #ifdef TESTING
    884             {
    885                 psImage *mask = psPixelsToMask(NULL, inspect->data[i],
    886                                                psRegionSet(0, numCols - 1, 0, numRows - 1),
    887                                                0xff); // Mask image
    888                 psString name = NULL;           // Name of image
    889                 psStringAppend(&name, "inspect_%03d.fits", i);
    890                 psFits *fits = psFitsOpen(name, "w");
    891                 psFree(name);
    892                 psFitsWriteImage(fits, NULL, mask, 0, NULL);
    893                 psFree(mask);
    894                 psFitsClose(fits);
    895             }
    896 #endif
    897 
    898             psPixels *reject = pmStackReject(inspect->data[i], numCols, numRows, threshold, poorFrac, stride,
    899                                              subRegions->data[i], subKernels->data[i]); // Rejected pixels
    900 
    901 #ifdef TESTING
    902             {
    903                 psImage *mask = psPixelsToMask(NULL, reject, psRegionSet(0, numCols - 1, 0, numRows - 1),
    904                                                0xff); // Mask image
    905                 psString name = NULL;           // Name of image
    906                 psStringAppend(&name, "reject_%03d.fits", i);
    907                 psFits *fits = psFitsOpen(name, "w");
    908                 psFree(name);
    909                 psFitsWriteImage(fits, NULL, mask, 0, NULL);
    910                 psFree(mask);
    911                 psFitsClose(fits);
    912             }
    913 #endif
    914 
    915             psFree(inspect->data[i]);
    916             inspect->data[i] = NULL;
    917 
    918             if (!reject) {
    919                 psWarning("Rejection on image %d didn't work --- reject entire image.", i);
    920                 numRejected++;
    921                 inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_REJECT;
    922             } else {
    923                 float frac = reject->n / (float)(numCols * numRows); // Pixel fraction
    924                 psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)",
    925                          reject->n, i, frac * 100.0);
    926                 if (frac > imageRej) {
    927                     psWarning("Image %d rejected completely because rejection fraction (%.3f) "
    928                               "exceeds limit (%.3f)", i, frac, imageRej);
    929                     psFree(reject);
    930                     // reject == NULL means reject image completely
    931                     reject = NULL;
    932                     inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_BAD;
    933                     numRejected++;
    934                 }
    935             }
    936 
    937             // Images without a list of rejected pixels (the list may be empty) are rejected completely
    938             rejected->data[i] = reject;
    939 
    940             if (stats) {
    941                 psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_REJECT", PS_META_DUPLICATE_OK,
    942                                  "Time to perform rejection", psTimerMark("PPSTACK_REJECT"));
    943                 psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_PIXELS", PS_META_DUPLICATE_OK,
    944                                  "Number of pixels rejected", reject ? reject->n : 0);
    945             }
    946             psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
    947                      psTimerClear("PPSTACK_REJECT"));
    948         }
    949 
    950         psFree(inspect);
    951         psFree(subKernels);
    952         psFree(subRegions);
    953 
    954         if (numRejected >= num - 1) {
    955             psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed.");
    956             psFree(stack);
    957             psFree(rejected);
    958             psFree(inputMask);
    959             psFree(view);
    960             psFree(outRO);
    961             psFree(covariances);
    962             psFree(matchChi2);
    963             return false;
    964         }
    965 
    966         if (stats) {
    967             psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_IMAGES", 0,
    968                              "Number of images rejected completely", numRejected);
    969         }
    970     }
    971 
    972     psLogMsg("ppStack", PS_LOG_INFO, "Stage 5 : Pixel Rejection : %f sec", psTimerClear("PPSTACK_STEPS"));
    973     psTimerStart("PPSTACK_FINAL");
    974 
    975     memDump("reject");
    97690
    97791    // Final combination
    97892    psTrace("ppStack", 2, "Final stack of convolved images....\n");
    979     {
    980         stack->lastScan = 0;            // Reset read
    981         bool status;                    // Status of read
    982         for (int numChunk = 0; true; numChunk++) {
    983             ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0);
    984             if (!status) {
    985                 // Something went wrong
    986                 psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
    987                 psFree(stack);
    988                 psFree(rejected);
    989                 psFree(inputMask);
    990                 psFree(view);
    991                 psFree(outRO);
    992                 psFree(covariances);
    993                 psFree(matchChi2);
    994                 return false;
    995             }
    996             if (!thread) {
    997                 // Nothing more to read
    998                 break;
    999             }
     93    if (!ppStackCombineFinal(stack, options, config)) {
     94        psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination.");
     95        psFree(stack);
     96        psFree(options);
     97        return false;
     98    }
     99    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
     100    ppStackMemDump("final");
    1000101
    1001             // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
    1002             psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
    1003             psArrayAdd(job->args, 1, thread);
    1004             psArrayAdd(job->args, 1, config);
    1005             psArrayAdd(job->args, 1, outRO);
    1006             psArrayAdd(job->args, 1, inputMask);
    1007             psArrayAdd(job->args, 1, rejected);
    1008             psArrayAdd(job->args, 1, weightings);
    1009             psArrayAdd(job->args, 1, matchChi2);
    1010             if (!psThreadJobAddPending(job)) {
    1011                 psFree(job);
    1012                 psFree(stack);
    1013                 psFree(rejected);
    1014                 psFree(inputMask);
    1015                 psFree(view);
    1016                 psFree(outRO);
    1017                 psFree(covariances);
    1018                 psFree(matchChi2);
    1019                 return false;
    1020             }
    1021             psFree(job);
    1022         }
    1023102
    1024         if (!psThreadPoolWait(true)) {
    1025             psError(PS_ERR_UNKNOWN, false, "Unable to do final combination.");
    1026             psFree(stack);
    1027             psFree(inputMask);
    1028             psFree(rejected);
    1029             psFree(view);
    1030             psFree(outRO);
    1031             psFree(covariances);
    1032             psFree(matchChi2);
    1033             return false;
    1034         }
     103    // Clean up
     104    psTrace("ppStack", 2, "Cleaning up after combination....\n");
     105    if (!ppStackCombineFinal(stack, options, config)) {
     106        psError(PS_ERR_UNKNOWN, false, "Unable to clean up.");
     107        psFree(stack);
     108        psFree(options);
     109        return false;
    1035110    }
     111    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
     112    ppStackMemDump("cleanup");
    1036113
    1037     memDump("final");
    1038 
    1039 #ifdef TESTING
    1040     writeImage("combined_final.fits", NULL, outRO->image, config);
    1041 #endif
    1042 
    1043     // Sum covariance matrices
    1044     for (int i = 0; i < num; i++) {
    1045         if (inputMask->data.U8[i]) {
    1046             psFree(covariances->data[i]);
    1047             covariances->data[i] = NULL;
    1048         }
    1049     }
    1050     outRO->covariance = psImageCovarianceSum(covariances);
    1051     psFree(covariances);
    1052     psFree(matchChi2);
    1053 
    1054     if (stats) {
    1055         psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_FINAL", 0,
    1056                          "Time to make final stack", psTimerMark("PPSTACK_FINAL"));
    1057     }
    1058 
    1059     psLogMsg("ppStack", PS_LOG_INFO, "Stage 6 : Final Stack : %f sec", psTimerClear("PPSTACK_STEPS"));
    1060     psLogMsg("ppStack", PS_LOG_INFO, "Time to make final stack: %f sec", psTimerClear("PPSTACK_FINAL"));
    1061 
    1062     psTrace("ppStack", 2, "Cleaning up after combination....\n");
    1063 
    1064 #if 0
    1065     // Ensure masked regions really look masked
    1066     {
    1067         psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
    1068         psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    1069         if (!pmReadoutMaskApply(outRO, maskBad)) {
    1070             psWarning("Unable to apply mask");
    1071         }
    1072     }
    1073 #endif
    1074 
    1075     // Close up
    1076     bool wcsDone = false;           // Have we done the WCS?
    1077     for (int i = 0; i < num; i++) {
    1078         if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    1079             continue;
    1080         }
    1081 
    1082         ppStackThread *thread = stack->threads->data[0]; // Representative stack
    1083         pmReadout *inRO = thread->readouts->data[i]; // Template readout
    1084         if (inRO && !wcsDone) {
    1085             // Copy astrometry over
    1086             wcsDone = true;
    1087             pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
    1088             pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
    1089             pmChip *outChip = outRO->parent->parent; // Output chip
    1090             pmFPA *outFPA = outChip->parent; // Output FPA
    1091             if (!outHDU || !inHDU) {
    1092                 psWarning("Unable to find HDU at FPA level to copy astrometry.");
    1093             } else {
    1094                 if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
    1095                     psErrorClear();
    1096                     psWarning("Unable to read WCS astrometry from input FPA.");
    1097                     wcsDone = false;
    1098                 } else {
    1099                     if (!outHDU->header) {
    1100                         outHDU->header = psMetadataAlloc();
    1101                     }
    1102                     if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
    1103                         psErrorClear();
    1104                         psWarning("Unable to write WCS astrometry to output FPA.");
    1105                         wcsDone = false;
    1106                     }
    1107                 }
    1108             }
    1109         }
    1110 
    1111         if (tempDelete) {
    1112             psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
    1113             psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
    1114             psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
    1115             if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
    1116                 unlink(varianceResolved) == -1) {
    1117                 psWarning("Unable to delete temporary files for image %d", i);
    1118             }
    1119             psFree(imageResolved);
    1120             psFree(maskResolved);
    1121             psFree(varianceResolved);
    1122         }
    1123     }
    1124     psFree(imageNames);
    1125     psFree(maskNames);
    1126     psFree(varianceNames);
    1127 
    1128     psFree(inputMask);
    1129114    psFree(stack);
    1130115
    1131     memDump("cleanup");
    1132116
    1133     // Generate binned JPEGs
    1134     {
    1135         int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level
    1136         int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level
    1137 
    1138         // Target cells
    1139         pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG1");
    1140         pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2");
    1141         psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    1142 
    1143         pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
    1144         if (!pmReadoutRebin(ro1, outRO, maskValue, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
    1145             psError(PS_ERR_UNKNOWN, false, "Unable to bin output.");
    1146             psFree(ro1);
    1147             psFree(ro2);
    1148             psFree(outRO);
    1149             return false;
    1150         }
    1151         psFree(ro1);
    1152         psFree(ro2);
    1153     }
    1154 
    1155     psLogMsg("ppStack", PS_LOG_INFO, "Stage 7 : WCS & JPEGS : %f sec", psTimerClear("PPSTACK_STEPS"));
    1156 
    1157     if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
    1158         psTrace("ppStack", 1, "Photometering stacked image....\n");
    1159 
    1160         psTimerStart("PPSTACK_PHOT");
    1161 
    1162         fileActivation(config, combineFiles, false);
    1163         fileActivation(config, photFiles, true);
    1164         pmFPAview *photView = filesIterateDown(config);
    1165         if (!ppStackPhotometry(config, outRO, photView)) {
    1166             psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
    1167             psFree(outRO);
    1168             psFree(photView);
    1169             return false;
    1170         }
    1171         psFree(photView);
    1172 
    1173         fileActivation(config, combineFiles, true);
    1174 
    1175         if (stats) {
    1176             pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File
    1177             pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
    1178             psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
    1179             psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", sources->n);
    1180             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0,
    1181                              "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
    1182         }
    1183         psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
    1184     }
    1185 
    1186     psLogMsg("ppStack", PS_LOG_INFO, "Stage 8 : Photometry Analysis : %f sec", psTimerClear("PPSTACK_STEPS"));
    1187 
    1188     psThreadPoolFinalize();
    1189 
    1190     memDump("phot");
    1191 
    1192     // Statistics on output
    1193     if (stats) {
    1194         psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
    1195         psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
    1196         psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    1197 
    1198         ppStatsFPA(stats, outRO->parent->parent->parent, view, maskBad, config);
    1199     }
    1200 
    1201     memDump("stats");
    1202 
    1203     // Put version information into the header
    1204     pmHDU *hdu = pmHDUFromCell(outRO->parent);
    1205     if (!hdu) {
    1206         psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
    1207         psFree(outRO);
    1208         psFree(view);
     117    // Photometry
     118    psTrace("ppStack", 1, "Photometering stacked image....\n");
     119    if (!ppStackPhotometry(options, config)) {
     120        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
     121        psFree(options);
    1209122        return false;
    1210123    }
    1211     if (!hdu->header) {
    1212         hdu->header = psMetadataAlloc();
    1213     }
    1214     ppStackVersionMetadata(hdu->header);
     124    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
     125    ppStackMemDump("photometry");
    1215126
    1216     psFree(outRO);
    1217     psFree(view);
    1218127
    1219     // Write out summary statistics
    1220     if (stats) {
    1221         psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_STACK", 0,
    1222                          "Time to do photometry", psTimerClear("PPSTACK_TOTAL"));
    1223 
    1224         const char *statsMDC = psMetadataConfigFormat(stats);
    1225         if (!statsMDC || strlen(statsMDC) == 0) {
    1226             psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
    1227         } else {
    1228             fprintf(statsFile, "%s", statsMDC);
    1229         }
    1230         psFree((void *)statsMDC);
    1231         fclose(statsFile);
    1232 
    1233         psFree(stats);
    1234     }
    1235 
    1236     // Write out the output files
    1237     if (!filesIterateUp(config)) {
     128    // Finish up
     129    psTrace("ppStack", 1, "Finishing up....\n");
     130    if (!ppStackFinish(options, config)) {
     131        psError(PS_ERR_UNKNOWN, false, "Unable to finish up.");
     132        psFree(options);
    1238133        return false;
    1239134    }
     135    psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
     136    ppStackMemDump("finish");
    1240137
    1241     psLogMsg("ppStack", PS_LOG_INFO, "Stage 9 : Final Output : %f sec", psTimerClear("PPSTACK_STEPS"));
    1242 
    1243     memDump("finish");
    1244 
     138    psFree(options);
    1245139    return true;
    1246140}
Note: See TracChangeset for help on using the changeset viewer.