IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19337


Ignore:
Timestamp:
Sep 3, 2008, 9:50:19 AM (18 years ago)
Author:
Paul Price
Message:

Threading stacks. Following the same scheme as for ppMerge, but cleaned up some. Rejection stage still needs to be threaded.

Location:
trunk/ppStack/src
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/Makefile.am

    r18346 r19337  
    1414        ppStackVersion.c        \
    1515        ppStackMatch.c          \
    16         ppStackSources.c
     16        ppStackSources.c        \
     17        ppStackThread.c
    1718
    1819noinst_HEADERS =                \

  • trunk/ppStack/src/ppStack.h

    r19214 r19337  
    55#define PPSTACK_INSPECT_PIXELS "PPSTACK.PIXELS" // Name of rejected pixels metadata items
    66
     7#include <pslib.h>
    78#include <psmodules.h>
     9
     10// Mask values for inputs
     11typedef enum {
     12    PPSTACK_MASK_MATCH  = 0x01,         // PSF-matching failed
     13    PPSTACK_MASK_REJECT = 0x02,         // Rejection failed
     14    PPSTACK_MASK_BAD    = 0x04,         // Bad image (too many pixels rejected)
     15    PPSTACK_MASK_ALL    = 0xff          // All errors
     16} ppStackMask;
     17
     18// Thread for stacking chunks
     19//
     20// Each input file contributes a readout, into which is read a chunk from that file
     21typedef struct {
     22    psArray *readouts;                  // Input readouts to read and stack
     23    bool read;                          // Has the scan been read?
     24    bool busy;                          // Is the scan being processed?
     25    int firstScan;                      // First row of the chunk to be read for this group
     26    int lastScan;                       // Last row of the chunk to be read for this group
     27} ppStackThread;
     28
     29// Allocator
     30ppStackThread *ppStackThreadAlloc(psArray *readouts // Inputs readouts to read and stack
     31    );
     32
     33// Data for threads
     34typedef struct {
     35    psArray *threads;                   // Threads doing stacking
     36    int lastScan;                       // Last row that's been read
     37    psArray *imageFits;                 // FITS file pointers for images
     38    psArray *maskFits;                  // FITS file pointers for masks
     39    psArray *weightFits;                // FITS file pointers for weights
     40} ppStackThreadData;
     41
     42// Set up thread data
     43ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, // Array of input cells
     44                                          const psArray *imageNames, // Names of images to read
     45                                          const psArray *maskNames, // Names of masks to read
     46                                          const psArray *weightNames, // Names of weight maps to read
     47                                          const pmConfig *config // Configuration
     48    );
     49
     50// Read chunk into the first available file thread
     51ppStackThread *ppStackThreadRead(bool *status, // Status of read
     52                                 ppStackThreadData *stack, // Stacks available for reading
     53                                 pmConfig *config, // Configuration
     54                                 int numChunk, // Chunk number (only for interest)
     55                                 int overlap // Overlap between subsequent scans
     56    );
     57
     58// Initialise the threads
     59void ppStackThreadInit(void);
    860
    961// Setup command-line arguments
    … …  
    3183
    3284// Perform stacking on a readout
    33 bool ppStackReadoutInitial(const pmConfig *config,   // Configuration
    34                     pmReadout *outRO,   // Output readout
    35                     const psArray *readouts, // Input readouts
    36                     const psArray *regions, // Array with array of regions used in each PSF matching
    37                     const psArray *kernels // Array with array of kernels used in each PSF matching
     85//
     86// Returns an array of pixels to inspect for each input image
     87psArray *ppStackReadoutInitial(const pmConfig *config,   // Configuration
     88                               pmReadout *outRO,   // Output readout
     89                               const psArray *readouts, // Input readouts
     90                               const psArray *regions, // Array with array of regions used in each PSF match
     91                               const psArray *kernels // Array with array of kernels used in each PSF match
     92    );
     93
     94// Thread entry point for ppStackReadoutInitial
     95bool ppStackReadoutInitialThread(psThreadJob *job // Job to process
    3896    );
    3997
    … …  
    43101                         const psArray *readouts, // Input readouts
    44102                         const psArray *rejected // Array with pixels rejected in each image
     103    );
     104
     105// Thread entry point for ppStackReadoutFinal
     106bool ppStackReadoutFinalThread(psThreadJob *job // Job to process
    45107    );
    46108
    … …  
    62124
    63125/// Convolve image to match specified seeing
    64 bool ppStackMatch(pmReadout *readout, ///< Readout to be convolved; replaced with output
     126bool ppStackMatch(pmReadout *readout, // Readout to be convolved; replaced with output
    65127                  psArray **regions, // Array of regions used in each PSF matching, returned
    66128                  psArray **kernels, // Array of kernels used in each PSF matching, returned
    67                   const psArray *sources, ///< Array of sources
    68                   const pmPSF *psf,     ///< Target PSF
    69                   psRandom *rng,        ///< Random number generator
    70                   const pmConfig *config ///< Configuration
     129                  const psArray *sources, // Array of sources
     130                  const pmPSF *psf,     // Target PSF
     131                  psRandom *rng,        // Random number generator
     132                  const pmConfig *config // Configuration
    71133    );
    72134
    … …  
    75137///
    76138/// Corrects the sources to have consistent magnitudes.  Returns the source lists.
    77 psArray *ppStackSourceListAdd(psArray *lists, ///< List to which to add, or NULL
    78                               psArray *sources, ///< Sources to add
    79                               const pmConfig *config ///< Configuration
     139psArray *ppStackSourceListAdd(psArray *lists, // List to which to add, or NULL
     140                              psArray *sources, // Sources to add
     141                              const pmConfig *config // Configuration
    80142    );
    81143
    … …  
    83145///
    84146/// Corrects the sources to have consistent magnitudes where possible.  Returns the sources
    85 psArray *ppStackSourceListCombine(psArray *lists, ///< Source lists
    86                                   const pmConfig *config ///< Configuration
     147psArray *ppStackSourceListCombine(psArray *lists, // Source lists
     148                                  const pmConfig *config // Configuration
    87149    );
    88150
    89151/// Print sources into a ds9 regions file
    90 void ppStackSourcesPrint(const psArray *sources ///< Sources to print
     152void ppStackSourcesPrint(const psArray *sources // Sources to print
    91153    );
    92154

  • trunk/ppStack/src/ppStackLoop.c

    r19283 r19337  
    1212#include "ppStack.h"
    1313
    14 #define TESTING
     14//#define TESTING
    1515
    1616#define WCS_TOLERANCE 0.001             // Tolerance for WCS
    … …  
    9595}
    9696
    97 #if 0
    98 // Set the data level for a list of files
    99 static void fileSetDataLevel(pmConfig *config, // Configuration
    100                              char **files, // Files for which to set level
    101                              pmFPALevel level // Level to set
    102                              )
    103 {
    104     assert(config);
    105     assert(files);
    106 
    107     for (int i = 0; files[i] != NULL; i++) {
    108         psArray *selected = pmFPAfileSelect(config->files, files[i]); // Selected files of interest
    109         for (int j = 0; j < selected->n; j++) {
    110             pmFPAfile *file = selected->data[j];
    111             assert(file);
    112             file->dataLevel = level;
    113         }
    114         psFree(selected);
    115     }
    116     return;
    117 }
    118 #endif
    119 
    12097// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
    12198static pmFPAview *filesIterateDown(pmConfig *config // Configuration
    … …  
    244221    }
    245222    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
    246     int numScans = psMetadataLookupS32(NULL, recipe, "ROWS"); // Number of scans to read at once
    247223
    248224    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    … …  
    412388    int numGood = 0;                    // Number of good frames
    413389    int numCols = 0, numRows = 0;       // Size of image
     390    psVector *inputMask = psVectorAlloc(num, PS_TYPE_U8); // Mask for inputs
     391    psVectorInit(inputMask, 0);
    414392    for (int i = 0; i < num; i++) {
    415393        psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
    … …  
    423401            psFree(targetPSF);
    424402            psFree(rng);
     403            psFree(inputMask);
    425404            return false;
    426405        }
    … …  
    439418            psFree(targetPSF);
    440419            psFree(rng);
     420            psFree(inputMask);
    441421            return false;
    442422        }
    … …  
    447427        if (!ppStackMatch(readout, &regions, &kernels, sources, targetPSF, rng, config)) {
    448428            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
     429            inputMask->data.U8[i] = PPSTACK_MASK_MATCH;
    449430            psErrorClear();
    450431            continue;
    … …  
    474455    if (numGood == 0) {
    475456        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
    476         return false;
    477     }
    478 
     457        psFree(subKernels);
     458        psFree(subRegions);
     459        psFree(cells);
     460        psFree(inputMask);
     461        return false;
     462    }
    479463
    480464#ifdef TESTING
    … …  
    482466#endif
    483467
     468    ppStackThreadInit();
     469    ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, weightNames,
     470                                                      config); // Data for stacking
     471    psFree(cells);
     472
    484473    memDump("preinitial");
    485474
    486     // Stack the convolved files
     475    // Stack the convolvd files
    487476    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
     477    pmReadout *outRO = NULL;            // Output readout
     478    pmFPAview *view = NULL;             // View to readout
    488479    {
     480        int row0, col0;                 // Offset for readout
     481        int numCols, numRows;           // Size of readout
     482        ppStackThread *thread = stack->threads->data[0]; // Representative thread
     483        if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
     484            psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
     485            psFree(subKernels);
     486            psFree(subRegions);
     487            psFree(stack);
     488            psFree(inputMask);
     489            return false;
     490        }
     491
    489492        pmFPAfileActivate(config->files, false, NULL);
    490493        fileActivation(config, combineFiles, true);
    491         pmFPAview *view = filesIterateDown(config);
     494        view = filesIterateDown(config);
    492495        if (!view) {
    493             psFree(cells);
    494496            psFree(subKernels);
    495497            psFree(subRegions);
    496             return false;
    497         }
     498            psFree(stack);
     499            psFree(inputMask);
     500            return false;
     501        }
     502
    498503        pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
    499         pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
    500 
    501         psArray *readouts = psArrayAlloc(num); // Readouts for convolved images
    502         for (int i = 0; i < num; i++) {
    503             pmCell *cell = cells->data[i]; // Cell of interest
    504             if (!cell) {
    505                 // Bad image
    506                 continue;
    507             }
    508             pmReadout *ro = cell->readouts->data[0]; // Readout of interest
    509             if (!ro) {
    510                 ro = pmReadoutAlloc(cell);
    511             }
    512             readouts->data[i] = ro; // Readout into which to read
    513         }
    514         psFree(cells);
    515 
    516         // FITS files
    517         psArray *imageFits  = psArrayAlloc(num);
    518         psArray *maskFits   = psArrayAlloc(num);
    519         psArray *weightFits = psArrayAlloc(num);
    520         for (int i = 0; i < num; i++) {
    521             if (!readouts->data[i]) {
    522                 // Bad image
    523                 continue;
    524             }
    525             // Resolved names
    526             psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
    527             psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
    528             psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
    529             imageFits->data[i] = psFitsOpen(imageResolved, "r");
    530             maskFits->data[i] = psFitsOpen(maskResolved, "r");
    531             weightFits->data[i] = psFitsOpen(weightResolved, "r");
    532             psFree(imageResolved);
    533             psFree(maskResolved);
    534             psFree(weightResolved);
    535             if (!imageFits->data[i] || !maskFits->data[i] || !weightFits->data[i]) {
    536                 psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s",
    537                         (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)weightNames->data[i]);
     504        outRO = pmReadoutAlloc(outCell); // Output readout
     505
     506        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     507        psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     508        if (!pmReadoutStackDefineOutput(outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
     509            psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
     510            psFree(subKernels);
     511            psFree(subRegions);
     512            psFree(stack);
     513            psFree(inputMask);
     514            psFree(view);
     515            psFree(outRO);
     516            return false;
     517        }
     518
     519        bool status;                    // Status of read
     520        for (int numChunk = 0; true; numChunk++) {
     521            ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap);
     522            if (!status) {
     523                // Something went wrong
     524                psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
    538525                psFree(subKernels);
    539526                psFree(subRegions);
     527                psFree(stack);
     528                psFree(inputMask);
     529                psFree(view);
     530                psFree(outRO);
    540531                return false;
    541532            }
    542         }
    543 
    544         // Read convolutions by chunks
    545         bool more = true;               // More to read?
    546         for (int numChunk = 0; more; numChunk++) {
    547             psTrace("ppStack", 2, "Initial stack of chunk %d....\n", numChunk);
     533            if (!thread) {
     534                // Nothing more to read
     535                break;
     536            }
     537
     538            // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
     539            psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
     540            psArrayAdd(job->args, 1, thread);
     541            psArrayAdd(job->args, 1, config);
     542            psArrayAdd(job->args, 1, outRO);
     543            psArrayAdd(job->args, 1, subRegions);
     544            psArrayAdd(job->args, 1, subKernels);
     545            if (!psThreadJobAddPending(job)) {
     546                psFree(job);
     547                psFree(subKernels);
     548                psFree(subRegions);
     549                psFree(stack);
     550                psFree(inputMask);
     551                psFree(view);
     552                psFree(outRO);
     553                return false;
     554            }
     555            psFree(job);
     556        }
     557
     558        if (!psThreadPoolWait(false)) {
     559            psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination.");
     560            psFree(subKernels);
     561            psFree(subRegions);
     562            psFree(stack);
     563            psFree(inputMask);
     564            psFree(view);
     565            psFree(outRO);
     566            return false;
     567        }
     568
     569        memDump("initial");
     570    }
     571
     572    // Pixel rejection
     573    psArray *rejected = psArrayAlloc(num);
     574    {
     575        psArray *inspect = psArrayAlloc(num); // Pixels to inspect
     576        int numRejected = 0;        // Number of inputs rejected completely
     577
     578        // Count images rejected out of hand
     579        for (int i = 0; i < num; i++) {
     580            if (inputMask->data.U8[i]) {
     581                numRejected++;
     582            }
     583            inspect->data[i] = psPixelsAllocEmpty(PIXELS_BUFFER);
     584        }
     585
     586        // Harvest the jobs, combining the inspection lists
     587        psThreadJob *job;           // Completed job
     588        while ((job = psThreadJobGetDone())) {
     589            psArray *results = job->results; // Results of job
     590            psAssert(results && results->n == num, "Results array.");
    548591            for (int i = 0; i < num; i++) {
    549                 pmReadout *readout = readouts->data[i];
    550                 if (!readout) {
    551                     // Bad image
     592                if (inputMask->data.U8[i]) {
    552593                    continue;
    553594                }
    554 
    555                 if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap,
    556                                         config) ||
    557                     !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap,
    558                                             config) ||
    559                     !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap,
    560                                               config)) {
    561                     psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
    562                     psFree(readouts);
    563                     psFree(subKernels);
    564                     psFree(subRegions);
    565                     psFree(outRO);
    566                     psFree(view);
    567                     return false;
    568                 }
    569             }
    570 
    571 #ifndef PS_NO_TRACE
    572             {
    573                 pmReadout *ro = NULL;   // Representative readout
    574                 for (int i = 0; i < num && !ro; i++) {
    575                     ro = readouts->data[i];
    576                 }
    577                 psAssert(ro, "There should be a readout here.");
    578                 psTrace("ppStack", 4, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
    579                         ro->row0, ro->row0 + ro->image->numRows);
    580             }
    581 #endif
    582 
    583             if (!ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)) {
    584                 psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
    585                 psFree(readouts);
    586                 psFree(subKernels);
    587                 psFree(subRegions);
    588                 psFree(outRO);
    589                 psFree(view);
    590                 return false;
    591             }
    592 
    593             for (int i = 0; i < num && more; i++) {
    594                 pmReadout *readout = readouts->data[i];
    595                 if (!readout) {
    596                     // Bad image
    597                     continue;
    598                 }
    599                 more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, config);
    600                 more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, config);
    601                 more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, config);
    602             }
    603         }
    604 
    605         memDump("initial");
    606 
    607         // Reset for the second read
    608         // Extract the rejection lists
    609         psArray *rejected = psArrayAlloc(num); // Pixels to inspect
    610         int numRejected = 0;            // Number of inputs rejected completely
     595                inspect->data[i] = psPixelsConcatenate(inspect->data[i], results->data[i]);
     596            }
     597            psFree(job);
     598        }
     599
    611600        for (int i = 0; i < num; i++) {
    612             pmReadout *ro = readouts->data[i]; // Readout of interest
    613             if (!ro) {
    614                 // Bad image
    615                 numRejected++;
     601            if (inputMask->data.U8[i]) {
    616602                continue;
    617603            }
    618             psPixels *inspect = psPixelsAllocEmpty(PIXELS_BUFFER); // Inspection list for this readout
    619             psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD,
    620                                                                "^" PPSTACK_INSPECT_PIXELS "$"); // Iterator
    621             psMetadataItem *item;
    622             while ((item = psMetadataGetAndIncrement(iter))) {
    623                 psPixels *pixels = item->data.V; // Rejected pixels
    624                 if (!pixels) {
    625                     continue;
    626                 }
    627                 psTrace("ppStack", 5, "Adding %ld pixels to inspect to image %d", pixels->n, i);
    628                 inspect = psPixelsConcatenate(inspect, pixels);
    629             }
    630             psFree(iter);
    631             psMetadataRemoveKey(ro->analysis, PPSTACK_INSPECT_PIXELS);
    632             pmReadoutFreeData(ro);
    633 
    634             psLogMsg("ppStack", PS_LOG_INFO, "%ld total pixels to inspect from image %d", inspect->n, i);
    635604
    636605#ifdef TESTING
    … …  
    648617#endif
    649618
    650             psPixels *reject = pmStackReject(inspect, numCols, numRows, threshold, poorFrac,
    651                                              subRegions->data[i], subKernels->data[i]); // Pixels to reject
     619            psPixels *reject = pmStackReject(inspect->data[i], numCols, numRows, threshold, poorFrac,
     620                                             subRegions->data[i], subKernels->data[i]); // Rejected pixels
    652621
    653622#ifdef TESTING
    … …  
    665634#endif
    666635
    667             psFree(inspect);
     636            psFree(inspect->data[i]);
     637            inspect->data[i] = NULL;
     638
    668639            if (!reject) {
    669640                psWarning("Rejection on image %d didn't work --- reject entire image.", i);
    670641                numRejected++;
     642                inputMask->data.U8[i] = PPSTACK_MASK_REJECT;
    671643            } else {
    672                 float frac = reject->n / (float)(outRO->image->numCols * outRO->image->numRows); // Pixel frac
     644                float frac = reject->n / (float)(numCols * numRows); // Pixel fraction
    673645                psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)",
    674                         reject->n, i, frac * 100.0);
     646                         reject->n, i, frac * 100.0);
    675647                if (frac > imageRej) {
    676648                    psWarning("Image %d rejected completely because rejection fraction (%.3f) "
    … …  
    679651                    // reject == NULL means reject image completely
    680652                    reject = NULL;
     653                    inputMask->data.U8[i] = PPSTACK_MASK_BAD;
    681654                    numRejected++;
    682655                }
    683656            }
    684657            rejected->data[i] = reject;
    685 
    686             memDump("reject");
    687 
    688         }
     658        }
     659
     660        psFree(inspect);
    689661        psFree(subKernels);
    690662        psFree(subRegions);
    … …  
    692664        if (numRejected == num) {
    693665            psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed.");
    694             psFree(readouts);
     666            psFree(stack);
    695667            psFree(rejected);
     668            psFree(inputMask);
     669            psFree(view);
    696670            psFree(outRO);
     671            return false;
     672        }
     673    }
     674
     675    memDump("reject");
     676
     677    // Final combination
     678    psTrace("ppStack", 2, "Final stack of convolved images....\n");
     679    {
     680        stack->lastScan = 0;            // Reset read
     681        bool status;                    // Status of read
     682        for (int numChunk = 0; true; numChunk++) {
     683            ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0);
     684            if (!status) {
     685                // Something went wrong
     686                psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
     687                psFree(stack);
     688                psFree(rejected);
     689                psFree(inputMask);
     690                psFree(view);
     691                psFree(outRO);
     692                return false;
     693            }
     694            if (!thread) {
     695                // Nothing more to read
     696                break;
     697            }
     698
     699            // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
     700            psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
     701            psArrayAdd(job->args, 1, thread);
     702            psArrayAdd(job->args, 1, config);
     703            psArrayAdd(job->args, 1, outRO);
     704            psArrayAdd(job->args, 1, rejected);
     705            if (!psThreadJobAddPending(job)) {
     706                psFree(job);
     707                psFree(stack);
     708                psFree(rejected);
     709                psFree(inputMask);
     710                psFree(view);
     711                psFree(outRO);
     712                return false;
     713            }
     714            psFree(job);
     715        }
     716
     717        if (!psThreadPoolWait(true)) {
     718            psError(PS_ERR_UNKNOWN, false, "Unable to do final combination.");
     719            psFree(stack);
     720            psFree(inputMask);
     721            psFree(rejected);
    697722            psFree(view);
    698             return false;
    699         }
    700 
    701         memDump("prefinal");
    702 
    703         // Read convolutions by chunks
    704         psTrace("ppStack", 2, "Final stack of convolved images....\n");
    705         more = true;
    706         for (int numChunk = 0; more; numChunk++) {
    707             psTrace("ppStack", 2, "Final stack of chunk %d....\n", numChunk);
    708             for (int i = 0; i < num; i++) {
    709                 if (!rejected->data[i]) {
    710                     continue;
    711                 }
    712                 pmReadout *readout = readouts->data[i];
    713                 assert(readout);
    714 
    715                 if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, 0, config) ||
    716                     !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, 0, config) ||
    717                     !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, 0,
    718                                               config)) {
    719                     psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
    720                     psFree(readouts);
    721                     psFree(rejected);
    722                     psFree(outRO);
    723                     psFree(view);
    724                     return false;
    725                 }
    726             }
    727 
    728 #ifndef PS_NO_TRACE
    729             {
    730                 pmReadout *ro = NULL;
    731                 for (int i = 0; i < num && !ro; i++) {
    732                     if (!rejected->data[i]) {
    733                         continue;
     723            psFree(outRO);
     724            return false;
     725        }
     726    }
     727
     728    memDump("final");
     729
     730    psTrace("ppStack", 2, "Cleaning up after combination....\n");
     731
     732    // Ensure masked regions really look masked
     733    {
     734        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
     735        psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     736        if (!pmReadoutMaskApply(outRO, maskBad)) {
     737            psWarning("Unable to apply mask");
     738        }
     739    }
     740
     741    // Close up
     742    bool wcsDone = false;           // Have we done the WCS?
     743    for (int i = 0; i < num; i++) {
     744        if (!inputMask->data.U8[i]) {
     745            continue;
     746        }
     747        if (!wcsDone) {
     748            // Copy astrometry over
     749            wcsDone = true;
     750            ppStackThread *thread = stack->threads->data[0]; // Representative stack
     751            pmReadout *inRO = thread->readouts->data[i]; // Template readout
     752            pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
     753            pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
     754            pmChip *outChip = outRO->parent->parent; // Output chip
     755            pmFPA *outFPA = outChip->parent; // Output FPA
     756            if (!outHDU || !inHDU) {
     757                psWarning("Unable to find HDU at FPA level to copy astrometry.");
     758            } else {
     759                if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
     760                    psErrorClear();
     761                    psWarning("Unable to read WCS astrometry from input FPA.");
     762                    wcsDone = false;
     763                } else {
     764                    if (!outHDU->header) {
     765                        outHDU->header = psMetadataAlloc();
    734766                    }
    735                     ro = readouts->data[i];
    736                 }
    737                 if (ro && ro->image) {
    738                     psTrace("ppStack", 4, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
    739                             ro->row0, ro->row0 + ro->image->numRows);
    740                 }
    741             }
    742 #endif
    743 
    744             if (!ppStackReadoutFinal(config, outRO, readouts, rejected)) {
    745                 psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
    746                 psFree(readouts);
    747                 psFree(rejected);
    748                 psFree(outRO);
    749                 psFree(view);
    750                 return false;
    751             }
    752 
    753             for (int i = 0; i < num && more; i++) {
    754                 pmReadout *readout = readouts->data[i];
    755                 if (!readout) {
    756                     continue;
    757                 }
    758                 more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, config);
    759                 more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, config);
    760                 more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, config);
    761             }
    762         }
    763 
    764         psFree(rejected);
    765         for (int i = 0; i < readouts->n; i++) {
    766             pmReadout *ro = readouts->data[i]; // Readout of interest
    767             if (!ro) {
    768                 continue;
    769             }
    770             pmReadoutFreeData(ro);
    771         }
    772         psFree(readouts);
    773 
    774         // Ensure masked regions really look masked
    775         {
    776             psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
    777             psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    778             if (!pmReadoutMaskApply(outRO, maskBad)) {
    779                 psWarning("Unable to apply mask");
    780             }
    781         }
    782 
    783         // Close up
    784         bool wcsDone = false;           // Have we done the WCS?
    785         for (int i = 0; i < num; i++) {
    786             if (!readouts->data[i]) {
    787                 continue;
    788             }
    789             if (!wcsDone) {
    790                 // Copy astrometry over
    791                 wcsDone = true;
    792                 pmReadout *inRO = readouts->data[i]; // Template readout
    793                 pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
    794                 pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
    795                 pmChip *outChip = outRO->parent->parent; // Output chip
    796                 pmFPA *outFPA = outChip->parent; // Output FPA
    797                 if (!outHDU || !inHDU) {
    798                     psWarning("Unable to find HDU at FPA level to copy astrometry.");
    799                 } else {
    800                     if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
     767                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
    801768                        psErrorClear();
    802                         psWarning("Unable to read WCS astrometry from input FPA.");
     769                        psWarning("Unable to write WCS astrometry to output FPA.");
    803770                        wcsDone = false;
    804                     } else {
    805                         if (!outHDU->header) {
    806                             outHDU->header = psMetadataAlloc();
    807                         }
    808                         if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
    809                             psErrorClear();
    810                             psWarning("Unable to write WCS astrometry to output FPA.");
    811                             wcsDone = false;
    812                         }
    813771                    }
    814772                }
    815773            }
    816 
    817             psFitsClose(imageFits->data[i]);
    818             psFitsClose(maskFits->data[i]);
    819             psFitsClose(weightFits->data[i]);
    820             imageFits->data[i] = NULL;
    821             maskFits->data[i] = NULL;
    822             weightFits->data[i] = NULL;
    823             if (tempDelete) {
    824                 psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
    825                 psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
    826                 psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
    827                 if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
    828                     unlink(weightResolved) == -1) {
    829                     psWarning("Unable to delete temporary files for image %d", i);
    830                 }
    831                 psFree(imageResolved);
    832                 psFree(maskResolved);
    833                 psFree(weightResolved);
    834             }
    835         }
    836         psFree(imageNames);
    837         psFree(maskNames);
    838         psFree(weightNames);
    839         psFree(imageFits);
    840         psFree(maskFits);
    841         psFree(weightFits);
    842 
    843         memDump("final");
    844 
    845         if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
    846 
    847             fileActivation(config, combineFiles, false);
    848             fileActivation(config, photFiles, true);
    849             pmFPAview *photView = filesIterateDown(config);
    850 
    851             psTrace("ppStack", 1, "Photometering stacked image....\n");
    852             if (!ppStackPhotometry(config, outRO, view)) {
    853                 psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
    854                 psFree(outRO);
    855                 psFree(photView);
    856                 psFree(view);
    857                 return false;
    858             }
     774        }
     775
     776        if (tempDelete) {
     777            psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
     778            psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
     779            psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
     780            if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
     781                unlink(weightResolved) == -1) {
     782                psWarning("Unable to delete temporary files for image %d", i);
     783            }
     784            psFree(imageResolved);
     785            psFree(maskResolved);
     786            psFree(weightResolved);
     787        }
     788    }
     789    psFree(imageNames);
     790    psFree(maskNames);
     791    psFree(weightNames);
     792
     793    psFree(inputMask);
     794    psFree(stack);
     795
     796    memDump("cleanup");
     797
     798    if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
     799        psTrace("ppStack", 1, "Photometering stacked image....\n");
     800
     801        fileActivation(config, combineFiles, false);
     802        fileActivation(config, photFiles, true);
     803        pmFPAview *photView = filesIterateDown(config);
     804        if (!ppStackPhotometry(config, outRO, photView)) {
     805            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
     806            psFree(outRO);
    859807            psFree(photView);
    860 
    861             fileActivation(config, combineFiles, true);
    862         }
    863 
    864         memDump("phot");
    865 
    866         // Statistics on output
    867         if (stats) {
    868             psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
    869             psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
    870             psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    871 
    872             ppStatsFPA(stats, outCell->parent->parent, view, maskBad, config);
    873         }
    874 
    875         memDump("stats");
    876 
    877         // Put version information into the header
    878         pmHDU *hdu = pmHDUFromCell(outRO->parent);
    879         if (!hdu) {
    880             psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
    881             return false;
    882         }
    883         if (!hdu->header) {
    884             hdu->header = psMetadataAlloc();
    885         }
    886         ppStackVersionMetadata(hdu->header);
    887 
    888         // Write out the output files
    889         filesIterateUp(config);
    890 
     808            return false;
     809        }
     810        psFree(photView);
     811
     812        fileActivation(config, combineFiles, true);
     813    }
     814
     815    memDump("phot");
     816
     817    // Statistics on output
     818    if (stats) {
     819        psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
     820        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
     821        psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     822
     823        ppStatsFPA(stats, outRO->parent->parent->parent, view, maskBad, config);
     824    }
     825
     826    memDump("stats");
     827
     828    // Put version information into the header
     829    pmHDU *hdu = pmHDUFromCell(outRO->parent);
     830    if (!hdu) {
     831        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
    891832        psFree(outRO);
    892833        psFree(view);
    893     }
     834        return false;
     835    }
     836    if (!hdu->header) {
     837        hdu->header = psMetadataAlloc();
     838    }
     839    ppStackVersionMetadata(hdu->header);
     840
     841    // Write out the output files
     842    filesIterateUp(config);
     843
     844    psFree(outRO);
     845    psFree(view);
    894846
    895847    // Write out summary statistics

  • trunk/ppStack/src/ppStackReadout.c

    r19267 r19337  
    1212//#define TESTING                  // Write debugging output?
    1313
    14 
    15 bool ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    16                            const psArray *regions, const psArray *kernels)
     14bool ppStackReadoutInitialThread(psThreadJob *job)
     15{
     16    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     17
     18    psArray *args = job->args;          // Arguments
     19    ppStackThread *thread = args->data[0]; // Thread
     20    pmConfig *config = args->data[1];   // Configuration
     21    pmReadout *outRO = args->data[2];   // Output readout
     22    psArray *subRegions = args->data[3]; // Regions for PSF-matching
     23    psArray *subKernels = args->data[4]; // Kernels for PSF-matching
     24
     25    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels);
     26
     27    job->results = inspect;
     28    thread->busy = false;
     29
     30    return inspect ? true : false;
     31}
     32
     33bool ppStackReadoutFinalThread(psThreadJob *job)
     34{
     35    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     36
     37    psArray *args = job->args;          // Arguments
     38    ppStackThread *thread = args->data[0]; // Thread
     39    pmConfig *config = args->data[1];   // Configuration
     40    pmReadout *outRO = args->data[2];   // Output readout
     41    psArray *rejected = args->data[3];  // Rejected pixels
     42
     43    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected); // Status of operation
     44
     45    thread->busy = false;
     46
     47    return status;
     48}
     49
     50
     51psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
     52                               const psArray *regions, const psArray *kernels)
    1753{
    1854    assert(config);
    … …  
    88124
    89125    // Save list of pixels to inspect
     126    psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
    90127    for (int i = 0; i < num; i++) {
    91128        pmStackData *data = stack->data[i]; // Data for this image
    … …  
    97134            continue;
    98135        }
    99         psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PPSTACK_INSPECT_PIXELS,
    100                          PS_DATA_PIXELS | PS_META_DUPLICATE_OK, "Pixels to inspect from initial combination",
    101                          data->inspect);
     136        inspect->data[i] = psMemIncrRefCounter(data->inspect);
    102137    }
    103138    psFree(stack);
    … …  
    105140    sectionNum++;
    106141
    107     return true;
     142    return inspect;
    108143}
    109144
Note: See TracChangeset for help on using the changeset viewer.