IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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.

File:
1 edited

Legend:

Unmodified
Added
Removed

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