IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 21, 2008, 1:33:36 PM (17 years ago)
Author:
bills
Message:

Added a program to replace streaks and to compare two image files
(comparsion not complete yet)
Split streaksremove.c up so to make it less wieldy and to facilitate
reusing the io functions.

Location:
trunk/magic/remove/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/magic/remove/src

    • Property svn:ignore

      •  

        old new  
        11streaksremove
         2streaksreplace
         3streakscompare

  • trunk/magic/remove/src/streaksremove.c

    r20810 r20816  
    11#include "streaksremove.h"
    2 #include "libgen.h"
    3 #include "unistd.h"
    4 
    5 static nebServer *ourNebServer = NULL;
    6 
    7 // to enhance clarity we don't propagate errors up the stack
    8 // we just bail out
    9 void streaksremoveExit(psString str, int exitCode) {
    10     psErrorStackPrint(stderr, str);
    11     exit(exitCode);
    12 }
    13 
    14 static void sFileFree(sFile *sfile)
    15 {
    16     psFree(sfile->resolved_name);
    17     psFree(sfile->name);
    18     psFree(sfile);
    19 }
    20 
    21 
    22 // getNebServer()
    23 //
    24 // it gets created the first time a nebulous file is accessed.
    25 // if config is passed in we are to create it and exit if not able.
    26 // if config is null we return NULL indicating that there are no nebulous files in use
    27 // The purpose of this is to avoid passing extra arguments between all of the functions.
    28 // but it is probably a bad idea.
    29 //
    30 static nebServer *getNebServer(pmConfig *config)
    31 {
    32     if (ourNebServer) {
    33         return ourNebServer;
    34     }
    35 
     2
     3static pmConfig *parseArguments(int argc, char **argv);
     4static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
     5static StreakPixels * getStreakPixels(streakFiles *sfiles, Streaks *streaks);
     6static void exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue);
     7static bool warpedPixel(streakFiles *sfiles, PixelPos *cellCoord);
     8static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue);
     9static void writeImages(streakFiles *sf, bool exciseImageCube);
     10static bool replicateOutputs(streakFiles *sfiles);
     11
     12int
     13main(int argc, char *argv[])
     14{
     15    bool status;
     16
     17    psLibInit(NULL);
     18
     19    pmConfig *config = parseArguments(argc, argv);
    3620    if (!config) {
    37         return NULL;
    38     }
    39 
    40     psString serverURI = getenv("NEB_SERVER");
    41 
    42     // XXX: Note: all of the flags to psError should be true, but I don't think nebclient
    43     // uses psError
    44     if (!serverURI) {
    45         bool status;
    46         serverURI = psMetadataLookupStr(&status, config->site, "NEB_SERVER");
    47         if (!status) {
    48             psError(PM_ERR_CONFIG, true, "failed to lookup config value for NEB_SERVER.");
    49             streaksremoveExit("", PS_EXIT_CONFIG_ERROR);
    50         }
    51     }
    52 
    53     if (!serverURI) {
    54         psError(PM_ERR_CONFIG, true, "Could not determine nebulous server URI.");
    55 
    56         streaksremoveExit("", PS_EXIT_CONFIG_ERROR);
    57     }
    58 
    59     ourNebServer = nebServerAlloc(serverURI);
    60     if (!ourNebServer) {
    61         psError(PM_ERR_SYS, true, "failed to create a nebServer object from %s.", serverURI);
    62         streaksremoveExit("", PS_EXIT_CONFIG_ERROR);
    63     }
    64 
    65     return ourNebServer;
    66 }
    67 
    68 static psString
    69 resolveFilename(pmConfig *config, sFile *sfile, bool create)
    70 {
    71     sfile->inNebulous = IN_NEBULOUS(sfile->name);
    72 
    73     if (sfile->inNebulous) {
    74         // make sure we have our neb server connection
    75         nebServer *server = getNebServer(config);
    76         if (create) {
    77             // delete the existing file, since there may be more than one
    78             // instance. It will get created below in pmConfigConvertFilename
    79             if (nebFind(server, sfile->name)) {
    80                 nebDelete(server, sfile->name);
     21        psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n");
     22        return PS_EXIT_CONFIG_ERROR;
     23    }
     24
     25    psMetadata *masks = psMetadataLookupMetadata(&status, config->recipes, "MASKS");
     26    if (!status) {
     27        psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n");
     28        return PS_EXIT_CONFIG_ERROR;
     29    }
     30    double maskStreak = (double) psMetadataLookupU8(&status, masks, "STREAK");
     31    if (!status) {
     32        psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n");
     33        return PS_EXIT_CONFIG_ERROR;
     34    }
     35
     36    psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
     37
     38    Streaks *streaks = readStreaksFile(streaksFileName);
     39    if (!streaks) {
     40        psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName);
     41        exit (PS_EXIT_PROG_ERROR);
     42    }
     43
     44    streakFiles *sfiles = openFiles(config, true);
     45    setupAstrometry(sfiles);
     46
     47    bool exciseAll = false;
     48    // --keepnonwarped is a test and debug mode
     49    bool keepNonWarpedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_WARPED");
     50
     51    // we need to check for non warped pixels unless we've been asked not to or the stage is diff
     52    // (By definition pixels in diff images were included in a diff)
     53    bool checkNonWarpedPixels = ! (keepNonWarpedPixels || (sfiles->stage == IPP_STAGE_DIFF) );
     54
     55    if (checkNonWarpedPixels ) {
     56        // From ICD:
     57        // In the raw and detrended images, the pixels which were not
     58        // included in any of the streak-processed warps must also be masked.
     59        // Note that the warp and difference skycells are only generated if more
     60        // than a small fraction of the pixels are lit by the input image.
     61        // if no skycells are provided sfiles->exciseAll is set to true
     62
     63        if (! computeWarpedPixels(sfiles) ) {
     64            // we have no choice to excise all pixels
     65            exciseAll = true;
     66        }
     67    }
     68   
     69    if (sfiles->stage == IPP_STAGE_RAW) {
     70        // copy PHU to output files
     71        copyPHU(sfiles, true);
     72
     73        // advance to the first image extension
     74        if (!streakFilesNextExtension(sfiles)) {
     75            psErrorStackPrint(stderr, "failed to advance to first extension of input images");
     76            exit (PS_EXIT_PROG_ERROR);
     77        }
     78    }
     79
     80    // Iterate through each component of the input (only one except for raw stage)
     81    do {
     82        bool exciseImageCube = false;
     83
     84        // read the images and copy the data from the inputs to the outputs
     85        // This also sets up the astrometry and initializes the recovery image(s)
     86
     87        if (!readAndCopyToOutput(sfiles, exciseAll)) {
     88            // false return value indicates that this extension is not an image and
     89            // the contents has already been copied to the output file.
     90            // we don't need to process this extesion for streaks.
     91            continue;
     92        }
     93
     94        if (!exciseAll) {
     95            // Identify pixels to mask because of streaks
     96
     97            StreakPixels *pixels = getStreakPixels(sfiles, streaks);
     98
     99            // XXX:
     100            //
     101            // PFS: Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
     102            // in a config somewhere?  Not sure how to properly set the maskValue.
     103
     104
     105            if (checkNonWarpedPixels) {
     106                exciseNonWarpedPixels(sfiles, maskStreak);
    81107            }
    82         }
    83     }
    84 
    85     return pmConfigConvertFilename(sfile->name, config, create, create);
    86 }
    87 
    88 sFile *sFileOpen(pmConfig *config, ippStage stage, psString fileSelect,
    89     psString outputFilename, bool required)
    90 {
    91     bool status;
    92     sFile *sfile = psAlloc(sizeof(sFile));
    93     memset(sfile, 0, sizeof(sFile));
    94 
    95     // We use psFits directly to read the image unless the stage is warp
    96     // or diff. In those cases we use the pmFPAfile functions to read the image file
    97     // to make managing the astrometry easy.
    98     // The reason we don't use pmFPAfile in all cases is that I was having trouble getting
    99     // all of the keywords in the raw image files written to the output destreaked files
    100 
    101     if (!CHIP_LEVEL_INPUT(stage) && !strcmp(fileSelect, "INPUT")) {
    102         // stage is warp or diff AND fileSelect eq "INPUT"
    103         // get data from pmFPAfile.
    104 
    105         // we need to know what the nebulous and real filenames are so we steal
    106         // some code from pmFPAfileDefineFromArgs
    107         // XXX: create a pmFPAfile function that does this and use it there
    108         psArray *infiles = psMetadataLookupPtr(&status, config->arguments, "INPUT");
    109         if (!status) {
    110             psError(PS_ERR_PROGRAMMING, false, "INPUT not found in config->arguments");
    111             streaksremoveExit("", PS_EXIT_PROG_ERROR);
    112         }
    113         if (infiles->n != 1) {
    114             psError(PS_ERR_IO, false, "Found n == %ld files in %s in arguments expencted 1\n",
    115                 infiles->n, "INPUT");
    116             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    117         }
    118         // end of file name lookup code adapted from pmFPAfileDefineFromArgs
    119         //
    120         sfile->name = psStringCopy(infiles->data[0]);
    121         sfile->inNebulous = IN_NEBULOUS(sfile->name);
    122 
    123         // XXX: I should probably be using a different file rule for diffs, but I don't
    124         // have an input rule that takes a diff image.
    125         // Since they're skycells, this should be compatible
    126         sfile->pmfile = pmFPAfileDefineFromArgs(&status, config, "PPSUB.INPUT", "INPUT");
    127         if (!sfile->pmfile) {
    128             psError(PS_ERR_IO, false, "Failed to define file for %s", fileSelect);
    129             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    130         }
    131 
    132         pmFPAview *view = pmFPAviewAlloc(0);
    133         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    134             psError(PS_ERR_UNKNOWN, false, "Failed to load input.");
    135             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    136         }
    137         psFree(view);
    138         sfile->resolved_name = psStringCopy(sfile->pmfile->filename);
    139        
    140         // copy header from fpu
    141         sfile->header = (psMetadata*) psMemIncrRefCounter(sfile->pmfile->fpa->hdu->header);
    142 
    143         return sfile;
    144     }
    145 
    146     // For all other files we use using psFits for i/o
    147 
    148     psString name = psMetadataLookupStr(&status, config->arguments, fileSelect);
    149     if (!status || !name) {
    150         if (required) {
    151             psError(PS_ERR_IO, false, "Failed to lookup name for %s", fileSelect);
    152             sFileFree(sfile);
    153             streaksremoveExit("", 1);
    154         }
    155         return NULL;
    156     }
    157 
    158     // if outputFilename is not null name it contains the "directory"
    159     // and outputFilename is the basename name of the file (or nebulous key)
    160     // and the file is to be opened for writing
    161     if (outputFilename) {
    162         psStringAppend(&sfile->name, "%s%s", name, outputFilename);
    163         sfile->resolved_name = resolveFilename(config, sfile, true);
    164         if (!sfile->resolved_name) {
    165             psError(PS_ERR_IO, false, "Failed to resolve filename for %s", sfile->name);
    166             sFileFree(sfile);
    167             streaksremoveExit("", 1);
    168         }
    169         sfile->fits = psFitsOpen(sfile->resolved_name, "w");
    170         if (sfile->fits) {
    171             sfile->fits->options = psFitsOptionsAlloc();
    172         }
    173     } else {
    174         sfile->name = psStringCopy(name);
    175         sfile->resolved_name = resolveFilename(config, sfile, false);
    176         if (!sfile->resolved_name) {
    177             psError(PS_ERR_IO, false, "Failed to resolve name for %s", sfile->name);
    178             sFileFree(sfile);
    179             streaksremoveExit("", 1);
    180         }
    181         sfile->fits = psFitsOpen(sfile->resolved_name, "r");
    182         if (sfile->fits) {
    183             sfile->nHDU = psFitsGetSize(sfile->fits);
    184         }
    185     }
    186 
    187     if (!sfile->fits) {
    188         psError(PS_ERR_IO, false, "failed to open fits file %s for %s",
    189                     sfile->resolved_name, outputFilename ? "writing" : "reading");
    190         sFileFree(sfile);
    191         streaksremoveExit("", 1);
    192     }
    193 
    194     return sfile;
    195 }
    196 
    197 
    198 static bool
    199 readAstrometry(streakFiles *sf)
    200 {
    201     pmHDU *phu = pmFPAviewThisPHU(sf->view, sf->inAstrom->fpa);
    202     if (phu) {
    203         bool status;
    204         char *ctype = psMetadataLookupStr(&status, phu->header, "CTYPE1");
    205         if (ctype) {
    206             sf->bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    207         }
    208     }
    209 
    210     if (sf->bilevelAstrometry) {
    211         // Do we get here for GPC1 ? I don't necessarily have an fpa for the input image
    212         if (!pmAstromReadBilevelMosaic(sf->inAstrom->fpa, phu->header)) {
    213             psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA.");
    214             return false;
    215         }
    216     } else {
    217         pmHDU *hdu = pmFPAviewThisHDU(sf->view, sf->inAstrom->fpa);
    218         PS_ASSERT_PTR_NON_NULL(hdu, 1)
    219         PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
    220 
    221         // we use a default FPA pixel scale of 1.0
    222         if (!pmAstromReadWCS (sf->inAstrom->fpa, sf->chip, hdu->header, 1.0)) {
    223             psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
    224             return false;
    225         }
    226     }
    227 
    228     return true;
    229 }
    230 
    231 // load the fpa containing the astrometry, find the appropriate chip and process the data
    232 static void
    233 setupAstromFromFPA(streakFiles *sf)
    234 {
    235     sf->view = pmFPAviewAlloc(0);
    236 
    237     if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) {
    238         psError(PS_ERR_UNKNOWN, false, "Failed to load input.");
    239         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    240     }
    241 
    242     while ((sf->chip = pmFPAviewNextChip(sf->view, sf->inAstrom->fpa, 1)))  {
    243         if (sf->inAstrom->fpa->chips->n == 1) {
    244             // There's only one chip in this FPA and we've got it so break
    245             break;
    246         }
    247         bool status;
    248         psString chip_name = psMetadataLookupStr(&status, sf->chip->concepts, "CHIP.NAME");
    249         if (!strcmp(chip_name, sf->class_id)) {
    250             break;
    251         }
    252     }
    253     if (!sf->chip) {
    254         psError(PS_ERR_UNKNOWN, true, "Failed to find chip with data.");
    255         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    256     }
    257     if (!pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_BEFORE)) {
    258         psError(PS_ERR_UNKNOWN, false, "failed to load chip");
    259         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    260     }
    261 
    262     readAstrometry(sf);
    263 }
    264 
    265 streakFiles *openFiles(pmConfig *config)
    266 {
    267     bool status;
    268     streakFiles *sf = psAlloc(sizeof(streakFiles));
    269     memset(sf, 0, sizeof(*sf));
    270 
    271     sf->config = config;
    272 
    273     // error checking is done by sFileOpen. If a file can't be opened we just exit
    274     ippStage stage = psMetadataLookupS32(&status, config->arguments, "STAGE");
    275 
    276     sf->stage = stage;
    277     sf->extnum = 0;
    278 
    279     sf->class_id = psMetadataLookupStr(&status, config->arguments, "CLASS_ID");
    280 
    281     sf->inImage = sFileOpen(config, stage, "INPUT", NULL, true);
    282     sf->nHDU = sf->inImage->nHDU;
    283 
    284     // don't need to free inputBasename see basename(3)
    285     // The names of the temporary and recovery files are taken from the input
    286     char *inputBasename = basename(sf->inImage->name);
    287     sf->outImage = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
    288     sf->recImage = sFileOpen(config, stage, "RECOVERY", inputBasename, false);
    289 
    290     sf->inMask = sFileOpen(config, stage, "INPUT.MASK", NULL, false);
    291     if (sf->inMask) {
    292         inputBasename = basename(sf->inMask->name);
    293         sf->outMask = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
    294         sf->recMask = sFileOpen(config, stage, "RECOVERY", inputBasename, true);
    295     }
    296     sf->inWeight = sFileOpen(config, stage, "INPUT.WEIGHT", NULL, false);
    297     if (sf->inWeight) {
    298         inputBasename = basename(sf->inWeight->name);
    299         sf->outWeight = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
    300         sf->recWeight = sFileOpen(config, stage, "RECOVERY", inputBasename, true);
    301     }
    302 
    303     // load astrometry file
    304     if (USE_SUPPLIED_ASTROM(sf->stage)) {
    305         sf->inAstrom = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM", "ASTROM");
    306     } else {
    307         // otherwise get astrometry from pmfile
    308         if (!sf->inImage->pmfile) {
    309             streaksremoveExit("unexpected null pmFPAfile", PS_EXIT_CONFIG_ERROR);
    310         }
    311         sf->inAstrom = sf->inImage->pmfile;
    312     }
    313     setupAstromFromFPA(sf);
    314      
    315     psElemType tileType;                // Type corresponding to "long"
    316     if (sizeof(long) == sizeof(psS64)) {
    317         tileType = PS_TYPE_S64;
    318     } else if (sizeof(long) == sizeof(psS32)) {
    319         tileType = PS_TYPE_S32;
    320     } else {
    321         psAbort("can't map (long) type to a psLib type");
    322     }
    323 
    324     sf->tiles = psVectorAlloc(3, tileType); // Tile sizes
    325     if (tileType == PS_TYPE_S64) {
    326         sf->tiles->data.S64[0] = 0;
    327         sf->tiles->data.S64[1] = 1;
    328         sf->tiles->data.S64[2] = 1;
    329     } else {
    330         sf->tiles->data.S32[0] = 0;
    331         sf->tiles->data.S32[1] = 1;
    332         sf->tiles->data.S32[2] = 1;
    333     }
    334 
    335     sf->transparentStreaks = psMetadataLookupF64(&status, config->arguments, "TRANSPARENT_STREAKS");
    336 
    337     return sf;
    338 }
    339 
    340 static bool
    341 moveExt(sFile *sfile, int extnum)
    342 {
    343     bool status = psFitsMoveExtNum(sfile->fits, extnum, false);
    344     if (!status) {
    345         psError(PS_ERR_IO, false,
    346             "failed to move to extension %d for %s", extnum, sfile->resolved_name);
    347         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    348     }
    349     return true;
    350 }
    351 
    352 static bool
    353 streakFilesNextExtension(streakFiles *sf)
    354 {
    355     // unless stage is raw or chip when we get here we're done
    356     if (sf->stage != IPP_STAGE_RAW) {
    357         sf->view->readout = -1;
    358         pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER);
    359         sf->view->cell = -1;
    360         pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER);
    361         sf->view->chip = -1;
    362         pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER);
    363         return false;
    364     }
    365 
    366     ++sf->extnum;
    367     if (sf->nHDU == 1) {
    368         // return true the first time through
    369         return (sf->extnum == 1);
    370     }
    371     if (sf->extnum < sf->nHDU) {
    372         moveExt(sf->inImage, sf->extnum);
    373         if (sf->inMask) {
    374             moveExt(sf->inMask, sf->extnum);
    375         }
    376         if (sf->inWeight) {
    377             moveExt(sf->inWeight, sf->extnum);
    378         }
    379         return true;
    380     } else {
    381         return false;
    382     }
    383 }
    384 
    385 psString checkdir(char *path)
    386 {
    387     // insure that path ends in /
    388     int len = strlen(path);
    389     psString dir = NULL;
    390     if (strrchr(path, '/') != (path + len - 1)) {
    391         psStringAppend(&dir, "%s/", path);
    392     } else {
    393         dir = psStringCopy(path);
    394     }
    395     return dir;
    396 }
    397 
    398 ippStage
    399 parseStage(psString stageStr)
    400 {
    401     if (!strcmp("raw", stageStr)) {
    402         return IPP_STAGE_RAW;
    403     } else if (!strcmp("chip", stageStr)) {
    404         return IPP_STAGE_CHIP;
    405     } else if (!strcmp("warp", stageStr)) {
    406         return IPP_STAGE_WARP;
    407     } else if (!strcmp("diff", stageStr)) {
    408         return IPP_STAGE_DIFF;
    409     } else {
    410         psError(PS_ERR_UNKNOWN, true, "unknown stage string: %s\n", stageStr);
    411         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    412         // notreached
    413         return IPP_STAGE_NONE;
    414     }
    415 }
    416        
    417 
    418 pmConfig *parseArguments(int argc, char **argv) {
     108           
     109            if (sfiles->inImage->image) {
     110                for (int i = 0; i < psArrayLength (pixels); ++i) {
     111                    PixelPos *pixelPos = psArrayGet (pixels, i);
     112
     113                    if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) {
     114
     115                        excisePixel(sfiles, pixelPos->x, pixelPos->y, true, maskStreak);
     116
     117                    } else {
     118                        // This pixel was not included in any warp and has thus already excised
     119                        // by exciseNonWarpedPixels
     120                    }
     121                }
     122            } else {
     123                // this component contains an image cube, excise it completely
     124                exciseImageCube = true;
     125            }
     126            psArrayElementsFree (pixels);
     127        }
     128
     129        // write out the destreaked temporary images and the recovery images
     130        writeImages(sfiles, exciseImageCube);
     131
     132    } while (streakFilesNextExtension(sfiles));
     133
     134    // close all files
     135    closeImages(sfiles);
     136
     137    // NOTE: from here on we can't just quit if something goes wrong.
     138    // especially if we're working at the raw stage
     139
     140    if (!replicateOutputs(sfiles)) {
     141        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
     142        psErrorStackPrint(stderr, "");
     143        exit(PS_EXIT_UNKNOWN_ERROR);
     144    }
     145
     146#ifdef NOTYET
     147    if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) {
     148        //     swap the instances for the input and output
     149        //     Note this is a database operation. No file I/O is performed
     150        if (!swapOutputsToInputs(sfiles)) {
     151            psError(PS_ERR_UNKNOWN, false, "failed to swap files");
     152
     153            // XXX: Now what? I guess swapOutputsToInputs will need to undo anything that
     154            // it has done and give a detailed report of what happened
     155
     156            psErrorStackPrint(stderr, "");
     157            exit(PS_EXIT_UNKNOWN_ERROR);
     158        }
     159
     160        if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) {
     161            //      delete the temporary storage objects (which now points to the original image(s)
     162            if (!deleteTemps(sfiles)) {
     163                psError(PS_ERR_UNKNOWN, false, "failed to delete temporary files");
     164                // XXX: Now what? At this point the output files have been swapped, so we can't
     165                // repeat the operation.
     166
     167                // Returning error status here is problematic. The inputs have been streak removed
     168                // but they're still lying around
     169                // Maybe just print an error message and
     170                // let other system tools clean up
     171                psErrorStackPrint(stderr, "");
     172                exit(PS_EXIT_UNKNOWN_ERROR);
     173            }
     174        }
     175    }
     176#endif  // REPLACE, REMOVE
     177    // nebServerFree(ourNebServer);
     178    psFree(config);
     179    pmConceptsDone();
     180    pmConfigDone();
     181    psLibFinalize();
     182
     183    //  PAU
     184    fprintf(stderr, "Found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "streaksremove");
     185
     186    return 0;
     187}
     188
     189static pmConfig *parseArguments(int argc, char **argv)
     190{
    419191
    420192    pmConfig *config = pmConfigRead(&argc, argv, NULL);
    421193    if (!config) {
    422         streaksremoveExit("failed to parse arguments", PS_EXIT_CONFIG_ERROR);
     194        streaksExit("failed to parse arguments", PS_EXIT_CONFIG_ERROR);
    423195    }
    424196
    … …  
    550322            return NULL;
    551323        }
    552         psString dir = checkdir(argv[argnum]);
     324        psString dir = pathToDirectory(argv[argnum]);
    553325        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "directory for temporary files",
    554326            dir);
    … …  
    561333    if ((argnum = psArgumentGet(argc, argv, "-recovery"))) {
    562334        psArgumentRemove(argnum, &argc, argv);
    563         psString dir = checkdir(argv[argnum]);
     335        psString dir = pathToDirectory(argv[argnum]);
    564336        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "directory for recovery files",
    565337            dir);
    … …  
    591363
    592364    return config;
    593 }
    594 
    595 static void
    596 copyPHU(streakFiles *sfiles)
    597 {
    598     psAssert(sfiles->stage == IPP_STAGE_RAW, "copyPHU should only be used for raw stage");
    599 
    600     psMetadata *imageHeader = psFitsReadHeader(NULL, sfiles->inImage->fits);
    601     if (!imageHeader) {
    602         psError(PS_ERR_IO, false, "failed to read primary header from %s\n",
    603             sfiles->inImage->resolved_name);
    604         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    605     }
    606 
    607     // TODO: add keyword indicating that streaks have been removed
    608     if (!psFitsWriteBlank(sfiles->outImage->fits, imageHeader, NULL)) {
    609         psError(PS_ERR_IO, false, "failed to write primary header to %s",
    610             sfiles->outImage->resolved_name);
    611         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    612     }
    613     // TODO: add keyword indicating that this is the recovery image
    614     if (sfiles->recImage && !psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) {
    615         psError(PS_ERR_IO, false, "failed to write primary header to %s",
    616             sfiles->recImage->resolved_name);
    617         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    618     }
    619     psFree(imageHeader);
    620 
    621     sFile *inMask = sfiles->inMask;
    622     if (inMask) {
    623         psMetadata *maskHeader = psFitsReadHeader(NULL, inMask->fits);
    624         if (!maskHeader) {
    625             psError(PS_ERR_IO, false, "failed to read primary header from %s\n",
    626                 sfiles->inMask->resolved_name);
    627             streaksremoveExit("", 1);
    628         }
    629         // TODO: add keyword indicating that streaks have been removed
    630         if (!psFitsWriteBlank(sfiles->outMask->fits, maskHeader, NULL)) {
    631             psError(PS_ERR_IO, false, "failed to write primary header to %s",
    632                 sfiles->outMask->resolved_name);
    633             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    634         }
    635         // TODO: add keyword indicating that this is the recovery image
    636         if (sfiles->recMask && !psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) {
    637             psError(PS_ERR_IO, false, "failed to write primary header to %s",
    638                 sfiles->recMask->resolved_name);
    639             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    640         }
    641         psFree(maskHeader);
    642     }
    643     sFile *inWeight = sfiles->inWeight;
    644     if (inWeight) {
    645         psMetadata *weightHeader = psFitsReadHeader(NULL, inWeight->fits);
    646         if (!weightHeader) {
    647             psError(PS_ERR_IO, false, "failed to read primary header from %s\n",
    648                 sfiles->inWeight->resolved_name);
    649             streaksremoveExit("", 1);
    650         }
    651         // TODO: add keyword indicating that streaks have been removed
    652         if (!psFitsWriteBlank(sfiles->outWeight->fits, weightHeader, NULL)) {
    653             psError(PS_ERR_IO, false, "failed to write primary header to %s",
    654                 sfiles->outWeight->resolved_name);
    655             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    656         }
    657         // TODO: add keyword indicating that this is a recovery image
    658         if (sfiles->recWeight && !psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) {
    659             psError(PS_ERR_IO, false, "failed to write primary header to %s",
    660                 sfiles->recWeight->resolved_name);
    661             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    662         }
    663         psFree(weightHeader);
    664     }
    665 }
    666 
    667 // Determine whether the current header for this file is an image or not
    668 // Find a cleaner way to do this
    669 static bool
    670 isImage(sFile *in)
    671 {
    672     bool status;
    673     psString exttype = psMetadataLookupStr(&status, in->header, "EXTTYPE");
    674     if (exttype && !strcmp(exttype, "IMAGE")) {
    675         return true;
    676     }
    677     // examine current header and determine if it is an image
    678     psString xtension = psMetadataLookupStr(&status, in->header, "XTENSION");
    679     if (xtension) {
    680         if (!strcmp(xtension, "IMAGE")) {
    681             return true;
    682         } else if (!strcmp(xtension, "BINTABLE")) {
    683             psString ztension = psMetadataLookupStr(&status, in->header, "ZTENSION");
    684             if (ztension && !strcmp(ztension, "IMAGE")) {
    685                 return true;
    686             }
    687         }
    688     } else if (psMetadataLookupBool(&status, in->header, "ZIMAGE")) {
    689         return true;
    690     } else if (in->nHDU == 1) {
    691         // no extensions in the file, can just return true? For now require
    692         // at least one dimension
    693         int naxis =  psMetadataLookupS32(&status, in->header, "NAXIS");
    694         return naxis > 0;
    695     }
    696     return false;
    697 }
    698 
    699 static void
    700 readImageFrom_pmFile(streakFiles *sf)
    701 {
    702     // XXX: This function assumes that it is only used for a single
    703     // chip single cell FPA (i.e. a skycell)
    704     pmFPAview *view = sf->view;
    705     assert(view->chip == 0);
    706     view->cell = 0;
    707     sf->cell =  pmFPAviewThisCell(view, sf->inImage->pmfile->fpa);
    708 
    709     if (!sf->cell) {
    710         streaksremoveExit("no cells found in chip", PS_EXIT_PROG_ERROR);
    711     }
    712     if (sf->cell->readouts->n != 1) {
    713         psError(PS_ERR_PROGRAMMING, true, "unexpected number of readouts found: %ld", sf->cell->readouts->n);
    714         streaksremoveExit("", PS_EXIT_PROG_ERROR);
    715     }
    716     view->readout = 0;
    717     pmReadout *readout = pmFPAviewThisReadout(view, sf->inImage->pmfile->fpa);
    718     if (!readout) {
    719         streaksremoveExit("readout 0 not found", PS_EXIT_PROG_ERROR);
    720     }
    721 
    722     sf->inImage->image = (psImage*) psMemIncrRefCounter(readout->image);
    723     sf->inImage->numCols = readout->image->numCols;
    724     sf->inImage->numRows = readout->image->numRows;
    725 }
    726 
    727 void
    728 setDataExtent(ippStage stage, sFile *in)
    729 {
    730     if (stage == IPP_STAGE_RAW) {
    731         psString datasec = psMetadataLookupStr(NULL, in->header, "DATASEC");
    732         if (!datasec) {
    733             psError(PS_ERR_IO, false, "failed to find DATASEC in header");
    734             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    735         }
    736         int xmin, xmax, ymin, ymax;
    737         sscanf(datasec, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax);
    738         in->numCols = xmax - xmin + 1;
    739         in->numRows = ymax - ymin + 1;
    740     } else {
    741         in->numCols = in->image->numCols;
    742         in->numRows = in->image->numRows;
    743     }
    744 }
    745 
    746 static void
    747 readImage(sFile *in, int extnum, ippStage stage)
    748 {
    749     psRegion region = {0, 0, 0, 0};
    750 
    751     if (in->header) psFree(in->header);
    752     if (in->image) {
    753         psFree(in->image);
    754         in->image = NULL;
    755     }
    756     if (in->imagecube) {
    757         psFree(in->imagecube);
    758         in->imagecube = NULL;
    759     }
    760 
    761     in->header = psFitsReadHeader(NULL, in->fits);
    762     if (!in->header) {
    763         psError(PS_ERR_IO, false, "failed to read header from %s extnum: %d",
    764             in->resolved_name, extnum);
    765         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    766     }
    767 
    768     if (!isImage(in)) {
    769         return;
    770     }
    771 
    772     in->numZPlanes = psMetadataLookupS32(NULL, in->header, "NAXIS3");
    773 
    774     if (in->numZPlanes == 0) {
    775         in->image = psFitsReadImage(in->fits, region, 0);
    776         if (!in->image) {
    777             psError(PS_ERR_IO, false, "failed to read image from %s extnum: %d",
    778                 in->resolved_name, extnum);
    779             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    780         }
    781     }  else {
    782         if (stage != IPP_STAGE_RAW) {
    783             psError(PS_ERR_PROGRAMMING, true, "unexpected image cube found in non-raw image");
    784             streaksremoveExit("", PS_EXIT_PROG_ERROR);
    785         }
    786         in->imagecube = psFitsReadImageCube(in->fits, region);
    787         if (!in->imagecube) {
    788             psError(PS_ERR_IO, false, "failed to read image cube from %s extnum: %d",
    789                 in->resolved_name, extnum);
    790             streaksremoveExit("", PS_EXIT_DATA_ERROR);
    791         }
    792         psImage *image = (psImage *) (in->imagecube->data[0]);
    793     }
    794     setDataExtent(stage, in);
    795 }
    796 
    797 static void
    798 setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero)
    799 {
    800     if (!sfile) {
    801         return;
    802     }
    803     if (sfile->fits->options) {
    804         psFree(sfile->fits->options);
    805     }
    806     sfile->fits->options = psFitsOptionsAlloc();
    807     sfile->fits->options->scaling = PS_FITS_SCALE_MANUAL;
    808     sfile->fits->options->bitpix = bitpix;
    809     sfile->fits->options->bscale = bscale;
    810     sfile->fits->options->bzero = bzero;
    811 }
    812 
    813 static void
    814 copyFitsOptions(sFile *out, sFile *rec, sFile *in)
    815 {
    816     // Get current BITPIX, BSCALE, BZERO, EXTNAME
    817     // Probably not necessary to look the numerical values up in this
    818     // way, but guards against changes to psLib and cfitsio FITS
    819     // handling.
    820     psMetadataItem *bitpixItem = psMetadataLookup(in->header, "BITPIX");
    821     psAssert(bitpixItem, "Every FITS image should have BITPIX");
    822     int bitpix = psMetadataItemParseS32(bitpixItem);
    823     psMetadataItem *bscaleItem = psMetadataLookup(in->header, "BSCALE");
    824 
    825     float bscale;
    826     if (!bscaleItem) {
    827         psWarning("BSCALE isn't set; defaulting to unity");
    828         bscale = 1.0;
    829     } else {
    830         bscale = psMetadataItemParseF32(bscaleItem);
    831     }
    832     psMetadataItem *bzeroItem = psMetadataLookup(in->header, "BZERO");
    833     float bzero;
    834     if (!bzeroItem) {
    835         psWarning("BZERO isn't set; defaulting to zero");
    836         bzero = 0.0;
    837     } else {
    838         bzero = psMetadataItemParseF32(bzeroItem);
    839     }
    840 
    841 #ifdef COMPRESS_OUTPUT
    842     setFitsOptions(out, bitpix, bscale, bzero);
    843     setFitsOptions(rec, bitpix, bscale, bzero);
    844 #endif
    845 }
    846 
    847 
    848 static void
    849 copyTable(sFile *out, sFile *in, int extnum)
    850 {
    851     bool mdok;
    852     psString xtension = psMetadataLookupStr(&mdok, in->header, "XTENSION");
    853     psString extname = psMetadataLookupStr(&mdok, in->header, "EXTNAME");
    854 
    855     if (!xtension) {
    856         psWarning("extnum %d has no image and xtension is not defined in %s",
    857             extnum, in->resolved_name);
    858         // XXX abort?
    859         return;
    860     }
    861     if (!extname) {
    862         psWarning("extnum %d has no image and extname not defined in %s",
    863             extnum, in->resolved_name);
    864         // XXX abort?
    865         return;
    866     }
    867     psArray *table = psFitsReadTable(in->fits);
    868     if (!table) {
    869         psError(PS_ERR_UNKNOWN, false, "failed to read table in extension %d from in->resolved name", extnum);
    870         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    871     }
    872 
    873     if (!psFitsWriteTable(out->fits, out->header, table, extname)) {
    874         psError(PS_ERR_UNKNOWN, false, "failed to copy table in extension %d", extnum);
    875         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    876     }
    877 }
    878 
    879 static void
    880 createNewImage(sFile *out, sFile *in, int extnum, double initValue)
    881 {
    882     out->image = psImageRecycle(out->image, in->image->numCols, in->image->numRows, in->image->type.type);
    883     if (!out->image) {
    884         psError(PS_ERR_UNKNOWN, false, "failed to allocate image for extnsion %d", extnum);
    885         streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR);
    886     }
    887     psImageInit(out->image, initValue);
    888 }
    889 
    890 static void
    891 setupImageRefs(sFile *out, sFile *recovery, sFile *in, int extnum, bool exciseAll)
    892 {
    893     if (!exciseAll) {
    894         // set output image to be the input image
    895         out->image = (psImage*) psMemIncrRefCounter(in->image);
    896         if (recovery) {
    897             // create a recovery image initialized to NAN
    898             createNewImage(recovery, in, extnum, NAN);
    899         }
    900     } else {
    901         // Excising all pixels in the input
    902         // set the output to an image all NANs
    903         createNewImage(out, in, extnum, NAN);
    904         if (recovery) {
    905             // save the whole input as the recovery image
    906             recovery->image = (psImage*) psMemIncrRefCounter(in->image);
    907         }
    908     }
    909365}
    910366
    … …  
    929385    if (!sf->astrom) {
    930386        psError(PS_ERR_UNKNOWN, false, "failed to set up astrometry for extension %d", sf->extnum);
    931         streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR);
     387        streaksExit("", PS_EXIT_UNKNOWN_ERROR);
    932388    }
    933389    if (sf->stage == IPP_STAGE_CHIP) {
    … …  
    936392        if (!pmAstromWriteWCS(sf->inImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
    937393            psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
    938             streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR);
     394            streaksExit("", PS_EXIT_UNKNOWN_ERROR);
    939395        }
    940396    }
    … …  
    958414
    959415    // set up the compression parameters
    960 #ifdef COMPRESS_OUTPUT
     416#ifdef STREAKS_COMPRESS_OUTPUT
    961417    // compression of the image pixels is disabled for now. Some consortium members
    962418    // have problems reading them
    … …  
    967423    // perhaps we should just use the definition of COMP_IMG in the configuration
    968424    psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
    969 #endif
    970425    if (sf->recImage) {
    971426        psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
    972427    }
     428#endif
    973429
    974430    if (sf->inMask) {
    … …  
    1010466    // area (no overscan)
    1011467
    1012 
    1013468    return true;
    1014469}
    1015470
    1016 static void
    1017 writeImage(sFile *sfile, psString extname, int extnum)
    1018 {
    1019     if (!sfile || !sfile->image) {
    1020         return;
    1021     }
    1022     if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) {
    1023         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
    1024             sfile->resolved_name, extnum);
    1025         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    1026     }
    1027     psFree(sfile->image);
    1028     sfile->image = NULL;
    1029     psFree(sfile->header);
    1030     sfile->header = NULL;
    1031 }
    1032 
    1033 static void
    1034 writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum)
    1035 {
    1036     if (!imagecube) {
    1037         return;
    1038     }
    1039     if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) {
    1040         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
    1041             sfile->resolved_name, extnum);
    1042         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    1043     }
    1044     psFree(sfile->header);
    1045     sfile->header = NULL;
    1046 }
     471
    1047472
    1048473static void
    … …  
    1098523}
    1099524
    1100 static void
    1101 closeImage(sFile *sfile)
    1102 {
    1103     if (!sfile) {
    1104         return;
    1105     }
    1106     if (sfile->fits && !psFitsClose(sfile->fits)) {
    1107         psError(PS_ERR_IO, false, "failed to close image to %s",
    1108             sfile->resolved_name);
    1109         streaksremoveExit("", PS_EXIT_DATA_ERROR);
    1110     }
    1111     psFree(sfile->header);
    1112     psFree(sfile->image);
    1113     psFree(sfile->imagecube);
    1114     psFree(sfile->name);
    1115     psFree(sfile->resolved_name);
    1116     psFree(sfile);
    1117 }
    1118 
    1119 static void
    1120 closeImages(streakFiles *sf)
    1121 {
    1122     closeImage(sf->inImage);
    1123     closeImage(sf->outImage);
    1124     closeImage(sf->recImage);
    1125     if (sf->inMask) {
    1126         closeImage(sf->inMask);
    1127         closeImage(sf->outMask);
    1128         closeImage(sf->recMask);
    1129     }
    1130     if (sf->inWeight) {
    1131         closeImage(sf->inWeight);
    1132         closeImage(sf->outWeight);
    1133         closeImage(sf->recWeight);
    1134     }
    1135 }
    1136 
    1137 static bool
    1138 replicate(sFile *sfile, void *xattr)
    1139 {
    1140     if (!sfile->inNebulous) {
    1141         return true;
    1142     }
    1143     nebServer *server = getNebServer(NULL);
    1144 
    1145     // for now just set "user.copies" to 2
    1146     if (!nebSetXattr(server, sfile->name, "user.copies", "2",  NEB_REPLACE)) {
    1147         psError(PM_ERR_UNKNOWN, true, "nebSetXattr failed for %s\n%s", sfile->name, nebErr(server));
    1148         return false;
    1149     }
    1150     if (!nebReplicate(server, sfile->name, NULL, NULL)) {
    1151         psError(PM_ERR_UNKNOWN, true, "nebSetXattr failed for %s\n%s", sfile->name, nebErr(server));
    1152         return false;
    1153     }
    1154     return true;
    1155 }
    1156 
    1157525// for any of the outputs that are stored in Nebulous set their extended attributes
    1158526// and replicate the files
    … …  
    1200568}
    1201569
    1202 static bool
    1203 swapOutputToInput(sFile *in, sFile *out)
    1204 {
    1205 
    1206     if (in->inNebulous) {
    1207         nebServer *server = getNebServer(NULL);
    1208 
    1209         if (!nebSwap(server, in->name, out->name)) {
    1210             psError(PM_ERR_SYS, true, "failed to swap files for: %s.",
    1211                 in->name, nebErr(server));
    1212             return false;
    1213         }
    1214     } else {
    1215         psString command = NULL;
    1216 
    1217         // regular files. use mv creating a backup
    1218         // NOTE: -b is a gnu option and may not always be available
    1219         psStringAppend(&command, "mv -b --suffix=.bak %s %s", out->resolved_name, in->resolved_name);
    1220         int status = system(command);
    1221         if (status != 0) {
    1222             psError(PM_ERR_SYS, true, "%s failed with status %d.",
    1223                 command, status);
    1224             psFree(command);
    1225             return false;
    1226         }
    1227         psFree(command);
    1228         // XXX: shouldn't I be doing mv in->resolved_name.bak out->resolved_name
    1229         // to complete the swap if !remove?
    1230     }
    1231 
    1232     return true;
    1233 }
    1234 static bool
    1235 swapOutputsToInputs(streakFiles *sfiles)
    1236 {
    1237     if (sfiles->inMask) {
    1238         if (!swapOutputToInput(sfiles->inMask, sfiles->outMask)) {
    1239             psError(PM_ERR_SYS, false, "failed to swap instances for Mask.");
    1240             return false;
    1241         }
    1242     }
    1243 
    1244     if (sfiles->inWeight) {
    1245         if (!swapOutputToInput(sfiles->inWeight, sfiles->outWeight)) {
    1246             psError(PM_ERR_SYS, false, "failed to swap instances for Weight.");
    1247             return false;
    1248         }
    1249     }
    1250 
    1251     if (!swapOutputToInput(sfiles->inImage, sfiles->outImage)) {
    1252         psError(PM_ERR_SYS, false, "failed to swap instances for Image.");
    1253         return false;
    1254     }
    1255 
    1256     return true;
    1257 }
    1258 
    1259 static bool
    1260 deleteFile(sFile *sfile)
    1261 {
    1262 
    1263     if (sfile->inNebulous) {
    1264         nebServer *server = getNebServer(NULL);
    1265         if (!nebDelete(server, sfile->name)) {
    1266             psError(PM_ERR_SYS, false, "failed to delete %s\n%s.", sfile->name,
    1267                 nebErr(server));
    1268             return false;
    1269         }
    1270     } else {
    1271         if (unlink(sfile->resolved_name)) {
    1272             psError(PM_ERR_SYS, false, "failed to delete %s\n%s.", sfile->resolved_name,
    1273                 strerror(errno));
    1274             return false;
    1275         }
    1276     }
    1277     return true;
    1278 }
    1279 
    1280 static bool
    1281 deleteTemps(streakFiles *sfiles)
    1282 {
    1283     if (sfiles->outMask) {
    1284         if (!deleteFile(sfiles->outMask)) {
    1285             psError(PM_ERR_SYS, false, "failed to delete Mask.");
    1286             return false;
    1287         }
    1288     }
    1289 
    1290     if (sfiles->outWeight) {
    1291         if (!deleteFile(sfiles->outWeight)) {
    1292             psError(PM_ERR_SYS, false, "failed to delete Weight.");
    1293             return false;
    1294         }
    1295     }
    1296 
    1297     if (!deleteFile(sfiles->outImage)) {
    1298         psError(PM_ERR_SYS, false, "failed to delete Image.");
    1299         return false;
    1300     }
    1301 
    1302     return true;
    1303 }
    1304570
    1305571static void
    … …  
    1312578
    1313579    double imageValue  = psImageGet (sfiles->inImage->image,  x, y);
    1314     if (sfiles->recImage) {
     580    if (sfiles->recImage && !isnan(imageValue) ) {
    1315581        psImageSet (sfiles->recImage->image,  x, y, imageValue);
    1316582    }
    … …  
    1327593    }
    1328594
    1329     // TODO:
    1330     // Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
    1331     // in a config somewhere?  Not sure how to properly set the maskValue.
    1332  
    1333595    if (sfiles->inWeight) {
    1334596        double weightValue = psImageGet (sfiles->inWeight->image, x, y);
    … …  
    1385647}
    1386648
    1387 bool
     649static bool
    1388650warpedPixel(streakFiles *sfiles, PixelPos *cellCoord)
    1389651{
    … …  
    1419681}
    1420682
    1421 bool
     683static bool
    1422684streakEndOnComponent(streak *streak, double r_min, double r_max, double d_min, double d_max)
    1423685{
    … …  
    1434696
    1435697
    1436 Streaks *
     698static Streaks *
    1437699restrictStreaksToComponent(streakFiles *sfiles, Streaks *streaks)
    1438700{
    … …  
    1494756    return pixels;
    1495757}
    1496 int
    1497 main(int argc, char *argv[])
    1498 {
    1499     long i;
    1500     bool status;
    1501     StreakPixels *pixels;
    1502     PixelPos *pixelPos;
    1503 
    1504     psLibInit(NULL);
    1505 
    1506     pmConfig *config = parseArguments(argc, argv);
    1507     if (!config) {
    1508         psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n");
    1509         return PS_EXIT_CONFIG_ERROR;
    1510     }
    1511 
    1512     psMetadata *masks = psMetadataLookupMetadata(&status, config->recipes, "MASKS");
    1513     if (!status) {
    1514         psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n");
    1515         return PS_EXIT_CONFIG_ERROR;
    1516     }
    1517     double maskStreak = (double) psMetadataLookupU8(&status, masks, "STREAK");
    1518     if (!status) {
    1519         psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n");
    1520         return PS_EXIT_CONFIG_ERROR;
    1521     }
    1522 
    1523     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
    1524 
    1525     Streaks *streaks = readStreaksFile(streaksFileName);
    1526     if (!streaks) {
    1527         psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName);
    1528         exit (PS_EXIT_PROG_ERROR);
    1529     }
    1530 
    1531     streakFiles *sfiles = openFiles(config);
    1532 
    1533     bool exciseAll = false;
    1534     // --keepnonwarped is a test and debug mode
    1535     bool keepNonWarpedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_WARPED");
    1536 
    1537     // we need to check for non warped pixels unless we've been asked not to or the stage is diff
    1538     // (By definition pixels in diff images were included in a diff)
    1539     bool checkNonWarpedPixels = ! (keepNonWarpedPixels || (sfiles->stage == IPP_STAGE_DIFF) );
    1540 
    1541     if (checkNonWarpedPixels ) {
    1542         // From ICD:
    1543         // In the raw and detrended images, the pixels which were not
    1544         // included in any of the streak-processed warps must also be masked.
    1545         // Note that the warp and difference skycells are only generated if more
    1546         // than a small fraction of the pixels are lit by the input image.
    1547         // if no skycells are provided sfiles->exciseAll is set to true
    1548 
    1549         if (! computeWarpedPixels(sfiles) ) {
    1550             // we have no choice to excise all pixels
    1551             exciseAll = true;
    1552         }
    1553     }
    1554    
    1555 
    1556     if (sfiles->stage == IPP_STAGE_RAW) {
    1557         // copy PHU to output files
    1558         copyPHU(sfiles);
    1559 
    1560         // advance to the first image extension
    1561         if (!streakFilesNextExtension(sfiles)) {
    1562             psErrorStackPrint(stderr, "failed to advance to first extension of input images");
    1563             exit (PS_EXIT_PROG_ERROR);
    1564         }
    1565     }
    1566 
    1567     // Iterate through components of input files
    1568     do {
    1569         bool exciseImageCube = false;
    1570 
    1571         // read the images and copy the data from the inputs to the outputs
    1572         // This also sets up the astrometry and initializes the recovery image(s)
    1573 
    1574         if (!readAndCopyToOutput(sfiles, exciseAll)) {
    1575             // false return value indicates that this extension is not an image and
    1576             // the contents has already been copied to the output file.
    1577             // we don't need to process this extesion for streaks.
    1578             continue;
    1579         }
    1580 
    1581         if (!exciseAll) {
    1582             // Identify pixels to mask because of streaks
    1583 
    1584             pixels = getStreakPixels(sfiles, streaks);
    1585 
    1586             // XXX:
    1587             //
    1588             // PFS: Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
    1589             // in a config somewhere?  Not sure how to properly set the maskValue.
    1590 
    1591 
    1592             if (checkNonWarpedPixels) {
    1593                 exciseNonWarpedPixels(sfiles, maskStreak);
    1594             }
    1595            
    1596             if (sfiles->inImage->image) {
    1597                 for (i = 0; i < psArrayLength (pixels); ++i) {
    1598                     pixelPos = psArrayGet (pixels, i);
    1599 
    1600                     if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) {
    1601 
    1602                         excisePixel(sfiles, pixelPos->x, pixelPos->y, true, maskStreak);
    1603 
    1604                     } else {
    1605                         // This pixel was not included in any warp and has thus already excised
    1606                         // by exciseNonWarpedPixels
    1607                     }
    1608                 }
    1609             } else {
    1610                 // this component contains an image cube, excise it completely
    1611                 exciseImageCube = true;
    1612             }
    1613             psArrayElementsFree (pixels);
    1614         }
    1615 
    1616         // write out the destreaked temporary images and the recovery images
    1617         writeImages(sfiles, exciseImageCube);
    1618 
    1619     } while (streakFilesNextExtension(sfiles));
    1620 
    1621     // close all files
    1622     closeImages(sfiles);
    1623 
    1624     // NOTE: from here on we can't just quit if something goes wrong.
    1625     // especially if we're working at the raw stage
    1626 
    1627     if (!replicateOutputs(sfiles)) {
    1628         psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
    1629         psErrorStackPrint(stderr, "");
    1630         exit(PS_EXIT_UNKNOWN_ERROR);
    1631     }
    1632 
    1633 #ifdef NOTYET
    1634     if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) {
    1635         //     swap the instances for the input and output
    1636         //     Note this is a database operation. No file I/O is performed
    1637         if (!swapOutputsToInputs(sfiles)) {
    1638             psError(PS_ERR_UNKNOWN, false, "failed to swap files");
    1639 
    1640             // XXX: Now what? I guess swapOutputsToInputs will need to undo anything that
    1641             // it has done and give a detailed report of what happened
    1642 
    1643             psErrorStackPrint(stderr, "");
    1644             exit(PS_EXIT_UNKNOWN_ERROR);
    1645         }
    1646 
    1647         if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) {
    1648             //      delete the temporary storage objects (which now points to the original image(s)
    1649             if (!deleteTemps(sfiles)) {
    1650                 psError(PS_ERR_UNKNOWN, false, "failed to delete temporary files");
    1651                 // XXX: Now what? At this point the output files have been swapped, so we can't
    1652                 // repeat the operation.
    1653 
    1654                 // Returning error status here is problematic. The inputs have been streak removed
    1655                 // but they're still lying around
    1656                 // Maybe just print an error message and
    1657                 // let other system tools clean up
    1658                 psErrorStackPrint(stderr, "");
    1659                 exit(PS_EXIT_UNKNOWN_ERROR);
    1660             }
    1661         }
    1662     }
    1663 #endif  // REPLACE, REMOVE
    1664     nebServer(ourNebServer);
    1665     psFree(config);
    1666     pmConceptsDone();
    1667     pmConfigDone();
    1668     psLibFinalize();
    1669 
    1670     //  PAU
    1671     fprintf(stderr, "Found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "streaksremove");
    1672 
    1673     return 0;
    1674 }
     758
Note: See TracChangeset for help on using the changeset viewer.