/*
 * streaksremove
 *
 * Convert satellite streak detctions into masks and remove the covered pixels from the
 * input images.
 * Optionally swap the inputs with the outputs so that subsequent references to the original
 * images access the destreaked versions.
 */

#include "streaksremove.h"

static pmConfig *parseArguments(int argc, char **argv);
static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
static long exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
static bool diffedPixel(streakFiles *sfiles, int x, int y);
static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue);
static void writeImages(streakFiles *sf, bool exciseImageCube);
static void updateAstrometry(streakFiles *sfiles);
static void censorSources(streakFiles *sfiles, psImageMaskType maskStreak);
static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonDiffedPixels, psU16 maskStreak);

char *streaksProgram = "streaksremove";

// Note: For clarity the flow of this program is in main().
// There is not a lot of error checking is done in main.
// Until the end, where we might be doing Nebulous operations, called functions exit when an error
// is encountered.
int
main(int argc, char *argv[])
{
    bool status;

    psLibInit(NULL);
    psTimerStart("TOTAL_TIME");

    pmConfig *config = parseArguments(argc, argv);
    if (!config) {
        psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n");
        streaksExit("", PS_EXIT_CONFIG_ERROR);
    }

    psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");

    // call Paul Sydney's code to parse the streaks file that DetectStreaks produced
    Streaks *streaks = readStreaksFile(streaksFileName);
    if (!streaks) {
        psError(PS_ERR_UNKNOWN, false, "failed to read streaks file: %s", streaksFileName);
        streaksExit("", PS_EXIT_PROG_ERROR);
    }

    // open all of the input and output files, save their descriptions in the streakFiles struct
    streakFiles *sfiles = openFiles(config, true, argv[0]);
    setupAstrometry(sfiles);
    sfiles->stats = psMetadataAlloc();

    // Optionally we can set pixels that are masked to NAN since they couldn't have been
    // examined for streaks. Usually this is done by the distribution system just prior
    // to release
    bool nanForRelease = psMetadataLookupBool(&status, config->arguments, "NAN_FOR_RELEASE");
    if (nanForRelease && (sfiles->inMask == NULL)) {
        psError(PS_ERR_UNKNOWN, true, "mask image must be supplied with -release.");
        streaksExit("", PS_EXIT_PROG_ERROR);
    }

    bool exciseAll = false;
    // --keepnondiffed is a test and debug mode
    bool keepNonDiffedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_DIFFED");

    // we need to check for non diffed pixels unless we've been asked not to or the stage is diff
    // (By definition pixels in diff images were included in the difference images)
    bool checkNonDiffedPixels = ! (keepNonDiffedPixels || (sfiles->stage == IPP_STAGE_DIFF) );

    if (checkNonDiffedPixels ) {
        // From magic ICD:
        // In the raw and detrended images, the pixels which were not
        // included in any of the streak-processed warps must also be masked.
        // Note that the warp and difference skycells are only generated if more
        // than a small fraction of the pixels are lit by the input image.
        // if no skycells are provided sfiles->exciseAll is set to true

        psTimerStart("COMPUTE_DIFFED_PIXELS");
        if (! computeDiffedPixels(sfiles) ) {
            // we have no choice to excise all pixels
            exciseAll = true;
        }
        psF64 cdp_t = psTimerClear("COMPUTE_DIFFED_PIXELS");
        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "COMPUTE_UNDIFFED_PIXELS", PS_META_REPLACE, "time to compute non-diffedpixels", cdp_t);
        psLogMsg("streaksremove", PS_LOG_INFO, "time to compute diffed pixels: %f\n", cdp_t);
    }

    if (sfiles->stage == IPP_STAGE_RAW) {
        // Except for raw stage, all of our (GPC1) files have one image extension.
        // Raw files have a phu and multiple extensions, one per chip
        // Since this is a raw file, copy it's PHU to output files
        copyPHU(sfiles, true, true);

        // advance to the first image extension
        if (!streakFilesNextExtension(sfiles)) {
            psErrorStackPrint(stderr, "failed to advance to first extension of input images");
            exit (PS_EXIT_PROG_ERROR);
        }
    }

    long totalPixels = 0;
    long totalStreakPixels = 0;
    long nonDiffedPixels = 0;

    // accumulators for the various timers
    psF64 gsp_t = 0;
    psF64 enw_t = 0;
    psF64 rms_t = 0;
    psF64 cs_t = 0;
    psF64 wi_t = 0;
    psF64 ua_t = 0;
    // Iterate through each component of the input (except for raw images there is only one)
    do {
        bool exciseImageCube = false;

        psTimerStart("PROCESS_COMPONENT");

        // read the images and copy the data from the inputs to the outputs
        // This also sets up the astrometry and initializes the recovery image(s)

        if (!readAndCopyToOutput(sfiles, exciseAll)) {
            // false return value indicates that this extension is not an image and
            // the contents has already been copied to the output file.
            // we don't need to process this extesion for streaks.
            continue;
        }

        // now that we've read the input files, lookup the mask values
        strkGetMaskValues(sfiles);

        totalPixels += sfiles->inImage->numRows * sfiles->inImage->numCols;

        if (!exciseAll) {
            // Identify pixels to mask because of streaks

            psTimerStart("GET_STREAK_PIXELS");

            // call Paul Sydney's code to compute the set of pixels that are covered by the detected streaks
            StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols,
                                                        sfiles->inImage->numRows);
            gsp_t +=  psTimerClear("GET_STREAK_PIXELS");
            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "GET_STREAK_PIXELS", PS_META_REPLACE, "", gsp_t);
            psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));

            // if this extension contained an image, excise the streaked pixels.
            // otherwise it contained an image cube (video cell) which is handled in the if block
            if (sfiles->inImage->image) {
                if (checkNonDiffedPixels) {
                    psTimerStart("EXCISE_NON_DIFFED");

                    // set non-diffed pixels and variance to NAN, mask to maskStreak (since the pixel
                    // is excised as part of the destreaking process)
                    nonDiffedPixels += exciseNonDiffedPixels(sfiles, sfiles->maskStreak);

                    enw_t +=  psTimerClear("EXCISE_NON_DIFFED");
                    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "EXCISE_NON_DIFFED", PS_META_REPLACE, "", enw_t);
                    psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non diffed pixels: %f\n", enw_t);
                }

                psTimerStart("REMOVE_STREAKS");

                totalStreakPixels += censorPixels(sfiles, pixels, checkNonDiffedPixels, sfiles->maskStreak);

                rms_t += psTimerClear("REMOVE_STREAKS");
                psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REMOVE_STREAKS", PS_META_REPLACE, "", enw_t);
                psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", rms_t);

                if (nanForRelease) {
                    // set any pixels that were masked, to NAN (unless they are already NAN)
                    setMaskedToNAN(sfiles, sfiles->maskMask, true);
                }

            } else {
                // this component contains an image cube
                // For now excise it completely
                exciseImageCube = true;
            }
            psFree(pixels);
        }

        if (sfiles->stage == IPP_STAGE_CHIP) {
            // as a convience to the user of the output, replace the bogus WCS transform in the
            // chip processed files with the data calcuated by psastro at the camera stage
            // (actually we use a linear approximation)
            psTimerStart("UPDATE_ASTROMETRY");
            updateAstrometry(sfiles);
            ua_t += psTimerClear("UPDATE_ASTROMETRY");
            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "UPDATE_ASTROMETRY", PS_META_REPLACE, "", ua_t);
            psLogMsg("streaksremove", PS_LOG_INFO, "time to update astrometry: %f\n", ua_t);
        }

        psTimerStart("CENSOR_SOURCES");
        censorSources(sfiles, sfiles->maskStreak);
        cs_t += psTimerClear("CENSOR_SOURCES");
        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CENSOR_SOURCES", PS_META_REPLACE, "", cs_t);

        // write the destreaked "temporary" images and the recovery images
        psTimerStart("WRITE_IMAGES");
        writeImages(sfiles, exciseImageCube);
        wi_t += psTimerClear("WRITE_IMAGES");
        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "WRITE_IMAGES", PS_META_REPLACE, "", wi_t);

        psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));

    } while (streakFilesNextExtension(sfiles));


    psFree(streaks);

    if (exciseAll) {
        totalStreakPixels = totalPixels;
    }

    psF64 streakFraction = (double) totalStreakPixels / totalPixels;
    psLogMsg("streaksremove", PS_LOG_INFO, "   total pixels:  %8ld\n", totalPixels);
    psLogMsg("streaksremove", PS_LOG_INFO, "  streak pixels:  %8ld %4.2f%%\n", totalStreakPixels, streakFraction * 100);
    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "STREAK_FRACTION", PS_META_REPLACE, "", streakFraction);

    psF64 nonDiffedFraction = (double) nonDiffedPixels / totalPixels;
    psLogMsg("streaksremove", PS_LOG_INFO, "nondiffed pixels:  %8ld %4.2f%%\n", nonDiffedPixels, nonDiffedFraction * 100);
    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "NONDIFFED_FRACTION", PS_META_REPLACE, "", nonDiffedFraction);

    // check the weight and mask files for extra extensions that might be in files
    // (covariance matrix for example)
    copyExtraExtensions(sfiles);

    // all done close the files. This is where the files are written so it can take a long time.

    psTimerStart("CLOSE_IMAGES");
    closeImages(sfiles);
    psF64 ci_t = psTimerClear("CLOSE_IMAGES");
    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CLOSE_IMAGES", PS_META_REPLACE, "", ci_t);

    psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", ci_t);

#ifdef DO_REPLICATE
    psTimerStart("REPLICATE_OUTPUTS");
    if (!replicateOutputs(sfiles)) {
        psErrorStackPrint(stderr, "failed to replicate output files");
        deleteTemps(sfiles);
        exit(PS_EXIT_UNKNOWN_ERROR);
    }
    psF64 ro_t = psTimerClear("REPLICATE_OUTPUTS");
    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REPLICATE_OUTPUTS", PS_META_REPLACE, "", ro_t);
#endif

    // NOTE: from here on we can't just quit if something goes wrong.
    // especially if we're working at the raw stage
    // turn off automatic deletion of output files by streaksExit
    setStreakFiles(NULL);

    if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) {
        //     swap the instances for the input and output
        //     Note this is a nebulous database operation. No file I/O is performed
        psTimerStart("SWAP_INSTANCES");
        if (!swapOutputsToInputs(sfiles)) {
            // It is up to the program that reverts failed destreak runs to insure that
            // any original non-destreaked input files that have been swapped are restored and that the de-streaked
            // versions are deleted.

            psErrorStackPrint(stderr, "failed to swap files");

            // XXX: pick a specific error code for this failure
            exit(PS_EXIT_UNKNOWN_ERROR);
        }
        psF64 si_t = psTimerClear("SWAP_INSTANCES");
        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "SWAP_INSTANCES", PS_META_REPLACE, "", si_t);
    }

    psF64 total_time = psTimerClear("TOTAL_TIME");
    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "TOTAL_TIME", PS_META_REPLACE, "", total_time);
    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", total_time);

    psString statsFileName= psMetadataLookupStr(&status, config->arguments, "STATS");

    if (statsFileName) {
        psString resolved = pmConfigConvertFilename(statsFileName, config, true, true); // Resolved filename
        if (!resolved) {
            psError(psErrorCodeLast(), false, "Unable to resolve statistics file name");
            return false;
        }
        if (!psMetadataConfigWrite(sfiles->stats, resolved, NULL)) {
            psError(psErrorCodeLast(), false, "Unable to serialize stats metadata.\n");
            psFree(resolved);
            return false;
        }
        psFree(resolved);
    }

    // all done. Clean up and look for memory leaks.

    psFree(sfiles);
    psFree(config);
    pmConceptsDone();
    pmModelClassCleanup();
    streaksNebulousCleanup();
    pmConfigDone();
    psLibFinalize();

    int leaks = psMemCheckLeaks (0, NULL, stdout, false);

    fprintf(stderr, "Found %d leaks at %s\n", leaks, "streaksremove");

    return 0;
}

static long
censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak)
{
    long streakPixels = 0;

    for (int y=0 ; y < sfiles->inImage->numRows; y++) {
        for (int x = 0; x < sfiles->inImage->numCols; x++) {
            if (psImageGet(pixels, x, y)) {
                if (!checkNonDiffedPixels || diffedPixel(sfiles, x, y)) {
                    ++streakPixels;

                    excisePixel(sfiles, x, y, true, maskStreak);

                } else {
                    // This pixel was not included in any warp and has thus already excised
                    // by exciseNonDiffedPixels
                }
            }
        }
    }
    return streakPixels;
}

static void usage()
{
    fprintf(stderr, "USAGE: streaksremove -stage raw|chip|warp|diff -image IMAGE.fits -tmproot directory -streaks streaks_file\n\n");
    fprintf(stderr, "Optional arguments:\n");
    fprintf(stderr, "\t-class_id class_id: class_id (required for raw and chip stages)\n");
    fprintf(stderr, "\t-recovery: directory to save excised pixel images\n");
    fprintf(stderr, "\t-mask MASK.fits mask file to de-streak\n");
    fprintf(stderr, "\t-weight WEIGHT.fits: weight file to de-streak\n");
    fprintf(stderr, "\t-replace: replace the input images with the output\n");
    fprintf(stderr, "\t-keepnondiffed: do not exise pixels that were not part of difference processing\n");
    fprintf(stderr, "\t-transparent val: instead of setting excicsed pixel to NAN add val\n");
    fprintf(stderr, "\t-chip_mask MASK.fits: name of mask for chip stage (camera stage mask is passed tih -mask)\n");

    exit(2);


}
static pmConfig *parseArguments(int argc, char **argv)
{

    pmConfig *config = pmConfigRead(&argc, argv, NULL);
    if (!config) {
        streaksExit("failed to parse arguments", PS_EXIT_CONFIG_ERROR);
    }

    int argnum;
    ippStage stage = IPP_STAGE_NONE;
    if ((argnum = psArgumentGet(argc, argv, "-stage"))) {
        psArgumentRemove(argnum, &argc, argv);
        stage = parseStage(argv[argnum]);
        psMetadataAddS32(config->arguments, PS_LIST_TAIL, "STAGE", 0,
                "pipeline stage for streak removal", stage);
        psArgumentRemove(argnum, &argc, argv);
    } else {
        psError(PS_ERR_UNKNOWN, true, "-stage is required\n");
        usage();
    }

    if ((argnum = psArgumentGet(argc, argv, "-class_id"))) {
        psArgumentRemove(argnum, &argc, argv);
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "CLASS_ID", 0,
                "class_id", argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    } else if ((stage == IPP_STAGE_RAW) || (stage == IPP_STAGE_CHIP)) {
        psError(PS_ERR_UNKNOWN, true, "-class_id is required\n");
        usage();
    }

    if ((argnum = psArgumentGet(argc, argv, "-streaks"))) {
        psArgumentRemove(argnum, &argc, argv);
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STREAKS", 0,
                "STREAKS", argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    } else {
        psError(PS_ERR_UNKNOWN, true, "-streaks is required\n");
        usage();
    }


    // If input image is in nebulous we require mask, weight, and tmproot to all be
    // nebulous paths
    bool nebulousImage = false;
    if ((argnum = psArgumentGet(argc, argv, "-image"))) {
        // for raw and chip level images we use psFits for I/O and ...
        nebulousImage = IN_NEBULOUS(argv[argnum+1]);
        if (CHIP_LEVEL_INPUT(stage)) {
            psArgumentRemove(argnum, &argc, argv);
            psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT", 0,
                    "name of input image", argv[argnum]);
            psArgumentRemove(argnum, &argc, argv);
        } else  {
            // ... for warped images we use pmFPAfiles
            if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "INPUT", "-image", NULL)) { ;
                psError(PS_ERR_UNKNOWN, false, "failed to process -image");
                usage();
            }
        }
    } else {
        psError(PS_ERR_UNKNOWN, true, "-image is required\n");
        usage();
    }

    bool gotReplace = false;
    if ((argnum = psArgumentGet(argc, argv, "-replace"))) {
        if (!nebulousImage) {
            psError(PS_ERR_UNKNOWN, true, "replace is only supported for nebulous files");
            usage();
        }
        gotReplace = true;
        psArgumentRemove(argnum, &argc, argv);
        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "REPLACE", 0, "replace input files",
            true);
    }

    if ((argnum = psArgumentGet(argc, argv, "-release"))) {
        psArgumentRemove(argnum, &argc, argv);
        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "NAN_FOR_RELEASE", 0, "set masked pixels to NAN",
            true);
    }

    if ((argnum = psArgumentGet(argc, argv, "-keepnondiffed"))) {
        psArgumentRemove(argnum, &argc, argv);
        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_DIFFED", 0,
            "skip excising of non diffed pixels", true);
    }

    if ((argnum = psArgumentGet(argc, argv, "-transparent"))) {
        if (gotReplace) {
            psError(PS_ERR_UNKNOWN, true, "cannot make transparent streaks if -replace");
            usage();
        }
        psArgumentRemove(argnum, &argc, argv);
        double transparentStreaks = atof(argv[argnum]);
        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0,
            "value to adjust excised pixels", transparentStreaks);
        psArgumentRemove(argnum, &argc, argv);
    }

    // if skycells are not provided then we have to execise all pixels  unless -keepnondiffed
    pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist");

    if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", NULL)) { ;
        if (CHIP_LEVEL_INPUT(stage)) {
            psError(PS_ERR_UNKNOWN, true, "-astrom is required for raw and chip stages\n");
            usage();
        }
    }

    if ((argnum = psArgumentGet(argc, argv, "-mask"))) {
        psArgumentRemove(argnum, &argc, argv);
        bool nebulousMask = IN_NEBULOUS(argv[argnum]);
        // This isn't a strict requirement if stage == IPP_STAGE_CHIP
        if (nebulousMask != nebulousImage) {
            psError(PS_ERR_UNKNOWN, true, "mask image must be %snebulous path with %s image path\n",
                nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular");
            usage();
        }
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.MASK", 0,
                "name of input mask image", argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    }

    if ((argnum = psArgumentGet(argc, argv, "-chip_mask"))) {
        if (stage != IPP_STAGE_CHIP) {
            psError(PS_ERR_UNKNOWN, true, "-chip_mask only allowed at chip stage");
            usage();
        }
        psArgumentRemove(argnum, &argc, argv);
        bool nebulousMask = IN_NEBULOUS(argv[argnum]);
        if (nebulousMask != nebulousImage) {
            psError(PS_ERR_UNKNOWN, true, "chip mask image must be %snebulous path with %s image path\n",
                nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular");
            usage();
        }
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.CHMASK", 0,
                "name of input chip mask image", argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    }

    if ((argnum = psArgumentGet(argc, argv, "-weight"))) {
        psArgumentRemove(argnum, &argc, argv);
        bool nebulousWeight = IN_NEBULOUS(argv[argnum]);
        if (nebulousWeight != nebulousImage) {
            psError(PS_ERR_UNKNOWN, true, "weight image must have %snebulous path with %s image path\n",
                nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular");
            usage();
        }
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.WEIGHT", 0,
                "name of input weight image", argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    }
    if ((argnum = psArgumentGet(argc, argv, "-sources"))) {
        psArgumentRemove(argnum, &argc, argv);
        bool nebulousSources = IN_NEBULOUS(argv[argnum]);
        if (nebulousSources != nebulousImage) {
            psError(PS_ERR_UNKNOWN, true, "sources file must have %snebulous path with %s image path\n",
                nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular");
            usage();
        }
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SOURCES", 0,
                "name of input sources file", argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    }

    if ((argnum = psArgumentGet(argc, argv, "-stats"))) {
        psArgumentRemove(argnum, &argc, argv);
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STATS", 0,
                "name of input stats file", argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    }

    if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) {
        psArgumentRemove(argnum, &argc, argv);
        bool nebulousTmp = IN_NEBULOUS(argv[argnum]);
        if (gotReplace && (nebulousTmp != nebulousImage)) {
            psError(PS_ERR_UNKNOWN, true, "tmproot must have %snebulous path with %s image path with -replace\n",
                nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular");
            usage();
        }
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "path for (temporary if replae) output files",
            argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    } else {
        psError(PS_ERR_UNKNOWN, true, "-tmproot is required\n");
        usage();
    }

    if ((argnum = psArgumentGet(argc, argv, "-recovery"))) {
        psArgumentRemove(argnum, &argc, argv);
        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "path for recovery files",
            argv[argnum]);
        psArgumentRemove(argnum, &argc, argv);
    } else if ((stage == IPP_STAGE_RAW) && gotReplace) {
        psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw with -replace\n");
        usage();
    }

    if (argc != 1) {
        psString unexpectedArguments = NULL;
        for (int i=1; i<argc; i++) {
            psStringAppend(&unexpectedArguments, "%s ", argv[i]);
        }
        psError(PS_ERR_UNKNOWN, true, "unexpected arguments: %s", unexpectedArguments);
        psFree(unexpectedArguments);
        usage();
    }

    return config;
}

static void
updateAstrometry(streakFiles *sf)
{
    if (sf->bilevelAstrometry) {
        if (!linearizeTransforms(sf->astrom)) {
            // fit failed, leave the transform in the file unchanged
            return;
        }
    }
    if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
        psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
        streaksExit("", PS_EXIT_UNKNOWN_ERROR);
    }
    if (sf->outMask) {
        pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
    }
    if (sf->outWeight) {
        pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
    }
}

static void
setStreakBits(psImage *maskImage, psU32 maskStreak)
{
    for (int y=0 ; y < maskImage->numRows; y++) {
        for (int x=0 ; x < maskImage->numCols; x++) {
            maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskStreak;
        }
    }
}

// set the image for the output files to the contents of the corresponding input file.
static bool
readAndCopyToOutput(streakFiles *sf, bool exciseAll)
{
    if (sf->inImage->pmfile) {
        // image data from pmFPAfile (diff or warp file)
        readImageFrom_pmFile(sf);
        sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa,
                sf->chip, true, sf->cell->concepts, sf->inImage->numCols, sf->inImage->numRows);
    } else {
        // image data directly from psFits
        readImage(sf->inImage, sf->extnum, sf->stage, false);
        if (SFILE_IS_IMAGE(sf->inImage)) {
            sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false,
                sf->inImage->header, sf->inImage->numCols, sf->inImage->numRows);
        }
    }
    if (!sf->astrom) {
        psError(PS_ERR_UNKNOWN, false, "failed to set up astrometry for extension %d", sf->extnum);
        streaksExit("", PS_EXIT_UNKNOWN_ERROR);
    }
    sf->outImage->header =  psMemIncrRefCounter(sf->inImage->header);
    if (sf->recImage) {
        sf->recImage->header = psMetadataCopy(NULL, sf->inImage->header);
        addRecoveryKeyword(sf->recImage->header);
    }
    addDestreakKeyword(sf->outImage->header, true);

    if (!SFILE_IS_IMAGE(sf->inImage)) {
        // fprintf(stderr, "skipping extnum: %d not an image\n", sf->extnum);
        // XXX: This extension is probably a video table.
        copyTable(sf->outImage, sf->inImage, sf->extnum);
        return false;
    }

    if (sf->inImage->image) {
        setupImageRefs(sf->outImage, sf->recImage, sf->inImage, sf->extnum, exciseAll);
    }  else {
        // Image cubes are handeled specially
    }

    // set up the compression parameters
    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage, sf->tiles);

    if (sf->inMask) {
        readImage(sf->inMask, sf->extnum, sf->stage, true);
        if (sf->outMask) {
            sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
            if (sf->recMask) {
                sf->recMask->header = psMetadataCopy(NULL, sf->outMask->header);
                addRecoveryKeyword(sf->recMask->header);
            }
            addDestreakKeyword(sf->outMask->header, true);
            // Note: we don't excise the mask pixels even if exciseAll is true.
            setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, false);
            if (exciseAll) {
                strkGetMaskValues(sf);

                // add the STREAK bit to the mask image pixels
                setStreakBits(sf->inMask->image, sf->maskStreak);
            }
            if (sf->outChMask) {
                sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header);
                sf->outChMask->image = (psImage *) psMemIncrRefCounter(sf->outMask->image);
                if (sf->recChMask) {
                    sf->recChMask->header = (psMetadata *) psMemIncrRefCounter(sf->recMask->header);
                    sf->recChMask->image = (psImage *) psMemIncrRefCounter(sf->recMask->image);
                }
            }

            copyFitsOptions(sf->outMask, sf->recMask, sf->inMask, sf->tiles);
            if (sf->outChMask) {
                copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask, sf->tiles);
            }
        }
    }

    if (sf->inWeight) {
        readImage(sf->inWeight, sf->extnum, sf->stage, false);
        sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
        if (sf->recWeight) {
            sf->recWeight->header = psMetadataCopy(NULL, sf->outWeight->header);
            addRecoveryKeyword(sf->recWeight->header);
        }
        addDestreakKeyword(sf->outWeight->header, true);
        setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll);

        copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight, sf->tiles);
    }

    return true;
}

static void
writeImages(streakFiles *sf, bool exciseImageCube)
{
    psString extname = NULL;
    if (sf->nHDU > 1) {
        bool mdok;
        extname = psMetadataLookupStr(&mdok, sf->inImage->header, "EXTNAME");
    }
    if (sf->inImage->numZPlanes == 0)  {
        // note exciseing complete images is handled in readAndCopyToOutput
        writeImage(sf->outImage, extname, sf->extnum);
        writeImage(sf->recImage, extname, sf->extnum);
    } else {
        // we have an image cube
        if (exciseImageCube) {
            // copy the entire input image to the recovery image
            writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
        } else {
            // otherwise write it to the output
            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
        }

        // Now deal with the other output image
        if (exciseImageCube) {
            // Set the values in the imagecube images to NAN and write them to the output image
            for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
                psImage *image = psArrayGet (sf->inImage->imagecube, i);
                // XXX: NAN isn't right. It should be the integer equivalent. Should use exciseValue
                // but it isn't set with this code path. Fix that.
                psImageInit(image, 65535);
            }
            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
        } else {
            if (sf->recImage) {
                // Set the values in the imagecube images to zero
                for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
                    psImage *image = psArrayGet (sf->inImage->imagecube, i);
                    psImageInit(image, 0);
                }
                // copy the entire zeroed image to the recovery image
                writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
            }
        }
        // XXX: check this
        // Assumption: there are no mask and weight images for video cells
        return;
    }
    if (sf->outMask) {
        writeImage(sf->outMask, extname, sf->extnum);
        writeImage(sf->recMask, extname, sf->extnum);
        if (sf->stage == IPP_STAGE_CHIP) {
            writeImage(sf->outChMask, extname, sf->extnum);
            writeImage(sf->recChMask, extname, sf->extnum);
        }
    }
    if (sf->outWeight) {
        writeImage(sf->outWeight, extname, sf->extnum);
        writeImage(sf->recWeight, extname, sf->extnum);
    }
}

static void
excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue)
{
    double exciseValue = sfiles->inImage->exciseValue;

    // we clip so that the streak calculation code doesn't have to
    if ((x < 0) || (x >= sfiles->inImage->numCols) || (y < 0) || (y >= sfiles->inImage->numRows)) {
        return;
    }

    if (sfiles->inImage->image->type.type == PS_TYPE_U16) {
        psU16 imageValue  = sfiles->inImage->image->data.U16[y][x];
        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
            sfiles->recImage->image->data.U16[y][x] = imageValue;
        }

        if (sfiles->transparentStreaks == 0) {
            sfiles->outImage->image->data.U16[y][x] = exciseValue;
        } else {
            if (streak) {
                // as a visualization aid don't mask the pixel, just change the intensity
                sfiles->outImage->image->data.U16[y][x] = imageValue + sfiles->transparentStreaks;
            } else {
                sfiles->outImage->image->data.U16[y][x] = exciseValue;
            }
        }
    } else if (sfiles->inImage->image->type.type == PS_TYPE_F32) {
        float imageValue  = sfiles->inImage->image->data.F32[y][x];
        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
            sfiles->recImage->image->data.F32[y][x] = imageValue;
        }

        if (sfiles->transparentStreaks == 0) {
            sfiles->outImage->image->data.F32[y][x] = exciseValue;
        } else {
            if (streak) {
                // as a visualization aid don't mask the pixel, just change the intensity
                sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
            } else {
                sfiles->outImage->image->data.F32[y][x] = exciseValue;
            }
        }
    } else {
        // We could handle F64 but I don't think we ever get them.
        // Should catch this earlier
        psError(PS_ERR_IO, true, "unexpected image type %x found", sfiles->inImage->image->type.type );
        streaksExit("", PS_EXIT_PROG_ERROR);
    }

    if (sfiles->outWeight) {
        // Assume that weight images are always a floating point type
        if (sfiles->inWeight->image->type.type == PS_TYPE_F32) {
            if (sfiles->recWeight) {
                sfiles->recWeight->image->data.F32[y][x] = sfiles->inWeight->image->data.F32[y][x];
            }
            sfiles->outWeight->image->data.F32[y][x] = NAN;
        } else if (sfiles->inWeight->image->type.type == PS_TYPE_F64) {
            if (sfiles->recWeight) {
                sfiles->recWeight->image->data.F64[y][x] = sfiles->inWeight->image->data.F64[y][x];
            }
            sfiles->outWeight->image->data.F64[y][x] = NAN;
        } else {
            // Should catch this earlier
            psError(PS_ERR_IO, true, "unexpected weight image type %x found", sfiles->inWeight->image->type.type );
            streaksExit("", PS_EXIT_PROG_ERROR);
        }
    }
    if (sfiles->outMask) {
        if (sfiles->recMask) {
            sfiles->recMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] =
                sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
        }
        sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= newMaskValue;
    }
}

static long
exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
{
    int cell_x0 = sfiles->astrom->cell_x0;
    int cell_y0 = sfiles->astrom->cell_y0;
    int xParity = sfiles->astrom->xParity;
    int yParity = sfiles->astrom->yParity;
    int numCols = sfiles->inImage->numCols; // for raw images this was calculated from the width of datasec
    int numRows = sfiles->inImage->numRows; // for raw images this was calculated from the height of datasec

    long excisedPixels = 0;

//    printf("%2d x0: %4d y0: %4d xpar: %d ypar: %d\n", sfiles->extnum, cell_x0, cell_y0, xParity, yParity);

    // XXX: in the case of raw images we are making an implicit assumption here
    // that the DATASEC begins at (x, y) = (0, 0) in the array for this image
    // While this is true for GPC1 it isn't necessarily true
    for (int yCell = 0; yCell < numRows; yCell++) {
        int yChip;
        if (yParity == 1) {
            yChip = cell_y0 + yCell;
        } else {
            yChip = cell_y0 - yCell;
        }

        psU8 *pixels = sfiles->diffedPixels->data.U8[yChip];

        if (xParity == 1) {
            pixels += cell_x0;
            for (int xCell = 0; xCell < numCols; xCell++, pixels++) {
                if (! *pixels ) {
                    excisePixel(sfiles, xCell, yCell, false, newMaskValue);
                    excisedPixels++;
                }
            }
        } else {
            // negative parity
            // point to the right most pixel in chip coordinates (lowest in cell coords)
            pixels += cell_x0 - numCols;
            for (int xCell = numCols - 1; xCell >= 0 ; xCell--, pixels++) {
                if (!*pixels) {
                    excisePixel(sfiles, xCell, yCell, false, newMaskValue);
                    excisedPixels++;
                }
            }
        }
    }
    return excisedPixels;
}

static bool
diffedPixel(streakFiles *sfiles, int x, int y)
{
    PixelPos chipCoord;

    if (!CHIP_LEVEL_INPUT(sfiles->stage)) {
        // if we're here on a skycell image by definition this pixel was diffed
        return true;
    }

    // we clip so that the streak calculation code doesn't have to
    // clipping here insures that we don't touch the overscan regions
    if ((x < 0) || (x >= sfiles->inImage->numCols) ||
        (y < 0) || (y >= sfiles->inImage->numRows)) {
        return false;
    }

    cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y);

    if (chipCoord.x < 0 || chipCoord.x >= sfiles->diffedPixels->numCols) {
        return false;
    }
    if (chipCoord.y < 0 || chipCoord.y >= sfiles->diffedPixels->numRows) {
        return false;
    }

    return psImageGet(sfiles->diffedPixels, chipCoord.x, chipCoord.y) ? true : false;
}

// read a sources file (.cmf) and remove any rows whose coordinate is convered by a
// streak mask
static void
censorSources(streakFiles *sfiles, psImageMaskType maskStreak)
{
    if ((!sfiles->inSources) || (!sfiles->outMask)) {
        return;
    }
    psImage *maskImage = sfiles->outMask->image;
    if (!maskImage) {
        psError(PS_ERR_IO, false, "maskImage is null!");
        streaksExit("", PS_EXIT_PROG_ERROR);
    }

    sFile *in = sfiles->inSources;
    sFile *out = sfiles->outSources;

    // Primary header, should be "something.hdr"
    {
        psMetadata *header = psFitsReadHeader(NULL, in->fits);
        if (!header) {
            psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }

        bool status;
        psString extname = psMetadataLookupStr(&status, header, "EXTNAME");
        if (!extname) {
            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }
        addDestreakKeyword(header, true);

        if (!psFitsWriteBlank(out->fits, header, extname)) {
            psError(PS_ERR_IO, false, "failed to write blank in header of %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }
        psFree(header);
    }

    // Extension with PSF fits, should be "CLASS_ID.psf"
    // This assumes that the sources extension comes first
    {
        if (!psFitsMoveExtNum(in->fits, 1, true)) {
            psErrorStackPrint(stderr, "failed to move to extension 1 in from %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }

        psMetadata *header = psFitsReadHeader(NULL, in->fits);
        if (!header) {
            psErrorStackPrint(stderr, "failed to read header from %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }
        psString extname = psMetadataLookupStr(NULL, header, "EXTNAME");
        if (!extname) {
            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }

#ifdef CHECK_MEM_USE
        bool printStats = false;
        size_t allocated = 0;
        size_t persistent = 0;
        size_t allocatedBefore = psMemStats(printStats, &allocated, &persistent);
        printf("Before: %12lld %12lld %12lld\n", (long long) allocated, (long long) persistent, (long long) allocatedBefore);
#endif

        psFitsTable *table = psFitsReadTableNew(in->fits);

#ifdef CHECK_MEM_USE
        size_t allocatedAfter = psMemStats(printStats, &allocated, &persistent);
        printf("After:  %12lld %12lld %12lld\n", (long long) allocated, (long long) persistent, (long long) allocatedAfter);
        printf("Change: %12lld\n", (long long) allocatedAfter - allocatedBefore);
#endif

        if (!table) {
            psErrorStackPrint(stderr, "failed to read table in %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }

        int numCensored = 0;
        bool *rowCensored = NULL;
        if (table->numRows > 0) {
            // For speed we get intimate with the table
            // avoid 2 * numRows lookups of the column number by looking
            // it up ahead of time and verifying that the type is what we
            // expect.
            int xCol = psFitsTableGetColumnNumber(table, "X_PSF");
            assert(xCol >= 0 && xCol < table->numCols);
            assert(table->columns[xCol].type == PS_DATA_F32);

            int yCol = psFitsTableGetColumnNumber(table, "Y_PSF");
            assert(yCol >= 0 && yCol < table->numCols);
            assert(table->columns[yCol].type == PS_DATA_F32);

            // We mark censored rows in this array.
            rowCensored = psAlloc(table->numRows * sizeof(bool));
            memset(rowCensored, 0, table->numRows * sizeof(bool));

            for (int i = 0 ; i < table->numRows; i++) {

#ifndef GO_SLOWER
                psF32 x = table->columns[xCol].data.F32[i];
                psF32 y = table->columns[yCol].data.F32[i];
#else
                psF32 x = psFitsTableGetF32(NULL, table, i, "X_PSF");
                psF32 y = psFitsTableGetF32(NULL, table, i, "Y_PSF");
#endif
                psImageMaskType mask;
                if ((x >= maskImage->numCols) || (y >= maskImage->numRows) || (x <  0) || (y < 0) || isnan(x) || isnan(y)) {
                    mask = maskStreak;
                } else {
                    mask = maskImage->data.PS_TYPE_IMAGE_MASK_DATA[(int)y][(int)x];
                }

                // Key the source if the center pixel is not masked with maskStreak
                if (mask & maskStreak) {
                    rowCensored[i] = true;
                    numCensored++;
                }
            }

            // remove the censored rows from the table
            if (!psFitsTableCensor(table, rowCensored)) {
                streaksExit("", PS_EXIT_PROG_ERROR);
            }

            addDestreakKeyword(header, true);
        }

        if (table->numCols > 0) {
            printf("Censored %d sources\n", numCensored);
            if (! psFitsWriteTableNew(out->fits, header, table, extname)) {
                psErrorStackPrint(stderr, "failed to write table to %s", out->resolved_name);
                streaksExit("", PS_EXIT_DATA_ERROR);
            }
        } else {
            printf("No input sources\n");
            // We'd like to write a blank table, but this is as good as it gets without more work to parse the
            // header and create dummy columns.
            if (!psFitsWriteBlank(out->fits, header, extname)) {
                psErrorStackPrint(stderr, "failed to write empty header to %s", out->resolved_name);
                streaksExit("", PS_EXIT_DATA_ERROR);
            }
        }
        psFree(header);
        psFree(table);
        psFree(rowCensored);
    }

    // Copy any other extensions in the file (deteff for example) to the output
    // XXX: This assumes that all we have is binary tables.
    for (int i = 2; i< in->nHDU; i++) {
        if (!psFitsMoveExtNum(in->fits, 1, true)) {
            psErrorStackPrint(stderr, "failed to move to extension number %d in %s", i, in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }
        psMetadata *header = psFitsReadHeader(NULL, in->fits);
        if (!header) {
            psErrorStackPrint(stderr, "failed to read header from %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }
        psString extname = psMetadataLookupStr(NULL, header, "EXTNAME");
        if (!extname) {
            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }

        psFitsTable *table = psFitsReadTableNew(in->fits);
        if (!table) {
            psErrorStackPrint(stderr, "failed to read table in %s", in->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }
        if (! psFitsWriteTableNew(out->fits, header, table, extname)) {
            psErrorStackPrint(stderr, "failed to write table to %s", out->resolved_name);
            streaksExit("", PS_EXIT_DATA_ERROR);
        }
        psFree(header);
        psFree(table);
    }

    if (!psFitsClose(out->fits)) {
        psErrorStackPrint(stderr, "failed to close table %s", out->resolved_name);
        streaksExit("", PS_EXIT_DATA_ERROR);
    }
    out->fits = NULL;   // freed by psFitsclose()
}
