#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <strings.h>

#include "ppstamp.h"
#include "pmAstrometryUtils.h"

#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS

typedef enum {
    PPSTAMP_OFF,
    PPSTAMP_PARTIALLY_ON,
    PPSTAMP_ON,
    PSTAMP_ERROR
} ppstampOverlap;

static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance);

// convert the input chip's transforms to the output
static bool convertWCS(pmHDU *outHDU, pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *roi)
{
    PS_ASSERT_PTR_NON_NULL(outHDU, 1);
    PS_ASSERT_PTR_NON_NULL(outHDU->header, 1);

    // Center of stamp in chip coordinates
    double center_x = (roi->x0 + roi->x1) / 2;
    double center_y = (roi->y0 + roi->y1) / 2;
    pmAstromObj *check = pmAstromObjAlloc();
    check->chip->x = center_x;
    check->chip->y = center_y;
    chipToSky(check, inFPA, inChip);
    double r_before = check->sky->r;
    double d_before = check->sky->d;
    psLogMsg("ppstampMakeStamp", 2, "Before fit  Center Pixel RA: %f   DEC: %f (degrees)\n",
	     r_before * PS_DEG_RAD, d_before * PS_DEG_RAD);

    if (outFPA->toSky->type != PS_PROJ_TAN) {
	// we need to make astrometry terms for the output which are centered on the output chip center,
	// but we keep the original plate scale
	psFree(outFPA->toSky);
	outFPA->toSky = psProjectionAlloc (r_before, d_before, inFPA->toSky->Xs, inFPA->toSky->Ys, PS_PROJ_TAN);

	if (!pmAstromLinearizeToSky(inFPA, inChip, outFPA, outChip, roi)) {
            psFree(check);
            psError(PS_ERR_UNKNOWN, false, "Failed to linearize astrometry\n");
            return false;
        }
        check->chip->x = (roi->x1 - roi->x0) / 2;
        check->chip->y = (roi->y1 - roi->y0) / 2;
        chipToSky(check, outFPA, outChip);
        double r_after = check->sky->r;
        double d_after = check->sky->d;

        psLogMsg("ppstampMakeStamp", 2, "After fit:  Center Pixel RA: %f   DEC: %f (degrees)\n",
		 r_after * PS_DEG_RAD, d_after * PS_DEG_RAD);
        psLogMsg("ppstampMakeStamp", 2, "Error in fit to astrometry %.2f    %.2f (arcseconds)\n", 
		 (r_after - r_before) *  PS_DEG_RAD * 3600, (d_after - d_before) * PS_DEG_RAD * 3600);

        // XXX: should we fail if the fit is bad ??

    } else {
        outChip->toFPA     = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0);
        outChip->fromFPA   = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50);
        if (!outChip->fromFPA) {
            psError(PS_ERR_UNKNOWN, false, "inversion of toFPA transform failed");
            return false;
        }
    }
    psFree(check);

    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to write WCS to output\n");
        return false;
    }

    return true;
}

static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options, pmAstromObj *center)
{
    pmChip    *outChip;
    pmFPAview *view = pmFPAviewAlloc(0);

    // our output file has a single chip
    view->chip = 0;
    outChip = pmFPAviewThisChip(view, output->fpa);
    psFree(view);

    // copy data from the input chip header to the output.
    // since some of the keywords might be duplicated we may not want to copy both

    pmHDU *inHDU  = pmHDUFromChip(inChip);
    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);

    if (inHDU->header) {
        outHDU->header = psMetadataCopy(outHDU->header, inHDU->header);
    } else if (!outHDU->header) {
        outHDU->header = psMetadataAlloc();
    }

    // copy the fpa concepts
    pmConceptsCopyFPA(output->fpa, input->fpa, false, false);
    pmConceptsCopyChip(outChip, inChip, false);
    pmCell *outCell = outChip->cells->data[0]; // The only output cell

    if (inChip->cells->n == 1) {
        pmConceptsCopyCell(outCell, inChip->cells->data[0]);
        // Need to fix up the trimsec and biassec to correspond to the output
        psMetadataItem *trimsec = psMetadataLookup(outCell->concepts, "CELL.TRIMSEC");
        psFree(trimsec->data.V);
        trimsec->data.V = NULL;
        psMetadataItem *biassec = psMetadataLookup(outCell->concepts, "CELL.BIASSEC");
        psFree(biassec->data.V);
        biassec->data.V = psListAlloc(NULL);
    } else {
        psList *inCells = psArrayToList(inChip->cells); // Input cells
        pmConceptsAverageCells(outCell, inCells, NULL, NULL, false);
        psFree(inCells);
    }

    // If input had WCS convert it for stamp
    if (input->fpa->toSky) {
        // copy the input fpa's transforms
        // XXX: if we change to making multiple stamps per process we'll need to copy
        // the transformations
        output->fpa->toTPA   = psMemIncrRefCounter(input->fpa->toTPA);
        output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA);
        output->fpa->toSky   = psMemIncrRefCounter(input->fpa->toSky);

        if (!convertWCS(outHDU, input->fpa, inChip, output->fpa, outChip, &options->roi)) {
            return false;
        }
    } else {
        psWarning("No WCS present in input\n");
    }

    psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "RA_DEG", PS_META_REPLACE, "Right Ascension of stamp center", RAD_TO_DEG(center->sky->r));
    psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "DEC_DEG", PS_META_REPLACE, "Declination of stamp center", RAD_TO_DEG(center->sky->d));

    ppstampVersionMetadata(outHDU->header, options);

    return true;
}

// Extract pixels from an image. If part of the bounds of the region are
// partially outside of the input those pixels are set to zero.
static psImage *extractStamp(psImage *image, psRegion region, double value)
{
    int width  = region.x1 - region.x0 + 0.5;
    int height = region.y1 - region.y0 + 0.5;

    if (width < 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative width\n");
        return NULL;
    }
    if (height < 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "negative height\n");
        return NULL;
    }

    int srcX  = region.x0;
    int srcY  = region.y0;
    int lastX = region.x1;
    int lastY = region.y1;
    if (lastX > (image->col0 + image->numCols)) {
        lastX = image->col0 + image->numCols;
    }
    if (lastY > (image->row0 + image->numRows)) {
        lastY = image->row0 + image->numRows;
    }

    int leftBlank = 0;
    int dstX  = 0;
    if (srcX < image->col0) {
        leftBlank = image->col0 - srcX;
        dstX += leftBlank;
        srcX = 0;
    }

    int copyWidth  = lastX - srcX;
    if (copyWidth <= 0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "copyWidth is invalid: %d\n", copyWidth);
        return NULL;
    }
    if (copyWidth > width) {
        copyWidth = width;
    }

    psImage *output = psImageAlloc(width, height, image->type.type);
    if (!output) {
        psError(PS_ERR_UNKNOWN, false, "psImageAlloc failed\n");
        return NULL;
    }

    if (!psImageInit(output, value)) {
        psError(PS_ERR_UNKNOWN, false, "psInitImage failed\n");
        return NULL;
    }

    psElemType  type = image->type.type;
    psS32       elementSize =  PSELEMTYPE_SIZEOF(type);
    int         copySize = copyWidth * elementSize;

    int     dstY = 0;
    psS32   skip = image->row0 - srcY;
    if (skip > 0) {
        dstY = skip;
        srcY = image->row0;
    }
    for ( ; dstY < height; srcY++, dstY++) {
        if (srcY >= lastY) {
            break;
        }
        // copy copyWidth pixels from srcX,srcY to dstX, dstY
        psU8 *pdst = output->data.U8[dstY] + (dstX * elementSize);

        psU8 *psrc = image->data.U8[srcY]  + (srcX * elementSize);

        memcpy(pdst, psrc, copySize);
    }

    return output;
}

// Build the postage stamp output file

static int makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input,
                pmChip *inChip, pmFPAview *view, pmAstromObj *center)
{
    int status = false;

    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTAMP.OUTPUT");
    if (!output) {
        psError(PS_ERR_UNKNOWN, false, "Can't find output data\n");
        return PS_EXIT_DATA_ERROR;
    }
    char *fpaName = psMetadataLookupStr(NULL, input->fpa->concepts, "FPA.OBS"); // Name of FPA
    pmFPAview *outview = pmFPAviewAlloc(0);
    pmFPAAddSourceFromView(output->fpa, fpaName, outview, output->format);

    outview->chip = 0;
    outview->cell = 0;
    pmCell *target =  pmFPAviewThisCell(outview, output->fpa); // Target cell
    psFree(outview);
    outview = NULL;

    //   psMetadataPrint(stderr, inChip->concepts, 0);

    // These default to zero. would that be ok?
    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE, "Binning in x", 1);
    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, "Binning in y", 1);

    pmFPAfile *srcFile;
    pmFPAview *srcView = pmFPAviewAlloc(0);
    if (psArrayLength(inChip->cells) > 1) {
        // we extract our postage stamp from a mosaic
        pmFPAfile *mosaic = ppstampBuildMosaic(config, input, view);
        if (mosaic == NULL) {
            return PS_EXIT_UNKNOWN_ERROR;
        }
        srcFile = mosaic;
        srcView->chip = 0;
    } else {
        srcFile = input;
        *srcView = *view;
    }

    // At this point we know we have one cell
    srcView->cell = 0;

    pmReadout *readout;
    while ((readout = pmFPAviewNextReadout (srcView, srcFile->fpa, 1)) != NULL) {
        if (!readout->data_exists) {
            psError(PS_ERR_UNKNOWN, false, "no data in input readout!\n");
            continue;
        }
        pmReadout *outReadout = pmReadoutAlloc(target);
        if (!outReadout) {
            psError(PS_ERR_UNKNOWN, false, "failed to allocate output readout\n");
            status = false;
            break;
        }

        psRegion extractRegion = options->roi;

        // Close your eyes while I hack around bug 986
        if (ppstampMegacamWorkaround) {
            // the coordinates of the mosaic are shifted 32 pixels from the chip
            // TODO does this always apply? For example I doubt that applies to
            // skycells.
            extractRegion.x0 -= 32;
            extractRegion.x1 -= 32;
        }

        outReadout->image = extractStamp(readout->image, extractRegion,  NAN);
        if (!outReadout->image) {
            psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp image\n");
            status = false;
            break;
        }
        if (readout->variance) {
            outReadout->variance = extractStamp(readout->variance, extractRegion,  NAN);
            if (!outReadout->variance) {
                psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp weight image\n");
                status = false;
                break;
            }
        }
        if (readout->mask) {
            outReadout->mask = extractStamp(readout->mask, extractRegion,  0);
            if (!outReadout->mask) {
                psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
                status = false;
                break;
            }

            if (!setMaskedToNAN(config, outReadout->image, outReadout->mask, outReadout->variance)) {
                 psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
                 status = false;
                 break;
            }
        }

        outReadout->data_exists = true;
        outReadout->parent->data_exists = true;
        outReadout->parent->parent->data_exists = true;
        status = true;

        psFree(outReadout); // drop reference
    }

    psFree(srcView);

    if (status) {
        status = copyMetadata(output, input, inChip, options, center);
    }
    return status ? PS_EXIT_SUCCESS : PS_EXIT_UNKNOWN_ERROR;
}



static bool regionContainsPoint(psRegion *r, psPlane *pt)
{
    if (pt->x < r->x0)
        return false;
    if (pt->x >= r->x1)
        return false;
    if (pt->y < r->y0)
        return false;
    if (pt->y >= r->y1)
        return false;

    return true;
}

// true if the inner region is equal to or completely contained in
// the outer region
static bool regionContainsRegion(psRegion *outer, psRegion *inner)
{
    if ((outer->x0 <= inner->x0) &&
        (outer->y0 <= inner->y0) &&
        (outer->x1 >= inner->x1) &&
        (outer->y1 >= inner->y1)) {
        return true;
    } else {
        return false;
    }
}


static void TPToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
{
    psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP);
    // convert from FP to chip
    psPlaneTransformApply(pt->chip, chip->fromFPA, pt->FP);
}
static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
{
    // convert from sky to TP to FP
    psProject(pt->TP, pt->sky, fpa->toSky);
    TPToChip(pt, fpa, chip);
}

static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
{
    // chip to FP
    psPlaneTransformApply(pt->FP, chip->toFPA, pt->chip);
    // FP to TP to sky
    psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP);
    psDeproject(pt->sky, pt->TP, fpa->toSky);
}

static void compareToBox(psRegion *region, psPlane *pt)
{
    if (pt->x < region->x0)
        region->x0 = pt->x;
    if (pt->x > region->x1)
        region->x1 = pt->x;
    if (pt->y < region->y0)
        region->y0 = pt->y;
    if (pt->y > region->y1)
        region->y1 = pt->y;
}

// For roi width and height given in sky coordinates find a bounding box in chip coordinates
// that encloses the requested range
static void findBoundingBox(ppstampOptions *options, pmFPA *fpa, pmChip *chip, pmAstromObj *center)
{
    pmAstromObj *pt = pmAstromObjAlloc();

    // calculate the four corners of the bounding box in sky coordinates, translate them to
    // chip coordinates and build the ROI by comparison.

    options->roi.x0 = INFINITY;
    options->roi.x1 = -INFINITY;
    options->roi.y0 = INFINITY;
    options->roi.y1 = -INFINITY;

    double dx = 0.5 * options->roip.dRA  / fpa->toSky->Xs;
    double dy = 0.5 * options->roip.dDEC / fpa->toSky->Ys;

    if (dx > 5000) {
        fprintf(stderr, "requested width %f too large reducing to 5000\n", dx);
        dx = 5000;
    }
    if (dy > 5000) {
        fprintf(stderr, "requested height %f too large reducing to 5000\n", dy);
        dy = 5000;
    }

    pt->TP->xErr = 0;
    pt->TP->yErr = 0;

    pt->TP->x = center->TP->x - dx;
    pt->TP->y = center->TP->y - dy;
    TPToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    pt->TP->x = center->TP->x + dx;
    pt->TP->y = center->TP->y - dy;
    TPToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    pt->TP->x = center->TP->x + dx;
    pt->TP->y = center->TP->y + dy;
    TPToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    pt->TP->x = center->TP->x - dx;
    pt->TP->y = center->TP->y + dy;
    TPToChip(pt, fpa, chip);
    compareToBox(&options->roi, pt->chip);

    psFree(pt);

    // save the width and height in case we need them later
    options->roip.dX = options->roi.x1 - options->roi.x0;
    options->roip.dY = options->roi.y1 - options->roi.y0;
}

// findROI
// calculate the region of interest in chip coordinates and determine whether that region
// is completely contained on a single chip

static ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view,
                    pmFPAfile *input, pmFPAfile *astrom, bool bilevelAstrometry,
                    pmChip *chip, pmAstromObj *center)
{
    psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
    psRegion    *chipBounds = ppstampChipRegion(chip);
    bool        onChip = false;
    ppstampOverlap   returnval = PPSTAMP_OFF;

//    fprintf(stderr, "ppstampChipBounds: %s\n", psRegionToString(*chipBounds));

    // set up the astrometry
    pmHDU *hdu = pmFPAviewThisHDU(view, astrom->fpa);
    PS_ASSERT_PTR_NON_NULL(hdu, 1)
    PS_ASSERT_PTR_NON_NULL(hdu->header, 1)

    // we can live without astrometry if the ROI is completely specified in pixel coordinates
    bool mustHaveAstrometry = false;
    if (options->roip.celestialCenter || options->roip.celestialRange) {
        mustHaveAstrometry = true;
    }
    if (bilevelAstrometry) {
        if (!pmAstromReadBilevelChip (chip, hdu->header) && mustHaveAstrometry) {
            psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel chip astrometry for input FPA.");
            return false;
        }
    } else {
        // we use a default FPA pixel scale of 1.0
        if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0) && mustHaveAstrometry) {
            psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
            return false;
        }
    }

    if (options->roip.celestialCenter) {

        center->sky->r = options->roip.centerRA;
        center->sky->d = options->roip.centerDEC;
        center->sky->rErr = 0;
        center->sky->dErr = 0;

        skyToChip(center, input->fpa, chip);

        if (regionContainsPoint(chipBounds, center->chip)) {
            psLogMsg("ppstampMakeStamp", 2, "Found center (%f %f) on chip: %s\n",
                center->chip->x, center->chip->y, chipName);
            onChip = true;
        }
    } else {
        // center specified in pixels.
        // If the user specified a name of a chip name wait until we get to that one.
        // If no chip name was specified, select this one (the first one that had data)
        if ((options->chipName == NULL) || !strcasecmp(chipName, options->chipName)) {
            psLogMsg("ppstampMakeStamp", 2, "Center on chip: %s\n", chipName);
            center->chip->x = options->roip.centerX;
            center->chip->y = options->roip.centerY;
            center->chip->xErr = 0;
            center->chip->yErr = 0;
            onChip = true;
            chipToSky(center, input->fpa, chip);
        }
    }

    if (onChip) {
        if (options->roip.celestialRange) {
            findBoundingBox(options, input->fpa, chip, center);
        } else {
            int width  = options->roip.dX;
            int height = options->roip.dY;
            if (width > 5000) {
                fprintf(stderr, "requested width %d too large reducing to 5000\n", width);
                width = 5000;
            }
            if (height > 5000) {
                fprintf(stderr, "requested height %d too large reducing to 5000\n", height);
                height = 5000;
            }

            // calculate the ROI in chip coordinates
            options->roi.x0 = center->chip->x - width / 2;
            options->roi.x1 = options->roi.x0 + width;
            options->roi.y0 = center->chip->y - height / 2;
            options->roi.y1 = options->roi.y0 + height;
        }


        if (regionContainsRegion(chipBounds, &options->roi)) {
            // returnval = findCell(view, options, input, center, ppCell);
            psLogMsg("ppstampMakeStamp", 2, "ROI contained on: %s\n", chipName);
            returnval = PPSTAMP_ON;
        } else {
            psLogMsg("ppstampMakeStamp", 2, "Partial Overlap chip %s\n", chipName);
            psLogMsg("ppstampMakeStamp", 2, "ROI:         %s\n", psRegionToString(options->roi));
            psLogMsg("ppstampMakeStamp", 2, "Chip Extent: %s\n", psRegionToString(*chipBounds));

            returnval = PPSTAMP_PARTIALLY_ON;
        }
    }
    psFree(chipBounds);

    return returnval;
}

int ppstampMakeStamp (pmConfig *config, ppstampOptions *options)
{
    bool        status = false;
    int        returnval = PS_EXIT_SUCCESS;;
    bool        foundOverlap = false;

    pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPSTAMP.INPUT");
    if (!status) {
        psError(PS_ERR_UNKNOWN, true, "Can't find input file!\n");
        return PS_EXIT_DATA_ERROR;
    }

    pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
    if (!astrom) {
        astrom = input;
    } else if (astrom->camera != input->camera) {
        psError(PS_ERR_UNKNOWN, true, "Input camera and astrometry camera do not match");
        return PS_EXIT_CONFIG_ERROR;
    }

    pmFPAview *view = pmFPAviewAlloc(0);// View for level of interest

    // files associated with the science image
    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to load input.");
        psFree (view);
        return PS_EXIT_DATA_ERROR;
    }
    bool bilevelAstrometry  = false;
    pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
    if (phu) {
        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
        if (ctype) {
            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
        }
    }
    if (bilevelAstrometry) {
        if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA.");
            psFree(view);
            return PS_EXIT_DATA_ERROR;
        }
    }

    // Loop over the chips and find the one that contains the center
    pmAstromObj *center = pmAstromObjAlloc();
    pmChip *chip;
    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
        bool allDone = false;

        if (!chip->process || !chip->file_exists) {
            continue;
        }

        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
            psError(PS_ERR_UNKNOWN, false, "failed to load chip");
            status = false;
            break;
        }

        ppstampOverlap overlap = findROI(options, view, input, astrom, bilevelAstrometry, chip, center);

        switch (overlap) {
        case PPSTAMP_OFF:
            // keep looking
            break;
        case PPSTAMP_ON:
        case PPSTAMP_PARTIALLY_ON:
            returnval = makeStamp(config, options, input, chip, view, center);
            allDone = true;
            foundOverlap = true;
            break;
        case PSTAMP_ERROR:
            returnval = PS_EXIT_UNKNOWN_ERROR;
            allDone = true;
            break;
        default:
            psError(PS_ERR_PROGRAMMING, false, "findROI returned unexpected value %d\n", overlap);
            break;
        }

        pmFPAfileIOChecks(config, view, PM_FPA_AFTER);

        if (allDone) {
            view->chip = -1;
            break;
        }
    }
    pmFPAfileIOChecks(config, view, PM_FPA_AFTER);

    psFree(center);
    psFree(view);

    if (!foundOverlap && (returnval == PS_EXIT_SUCCESS)) {
        fprintf(stderr, "ROI not found in input\n");
        returnval = PSTAMP_NO_OVERLAP;
    }

    return returnval;
}


static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance)
{
    bool status;
    psMetadata *masks = psMetadataLookupMetadata(&status, config->recipes, "MASKS");
    if (!status) {
        psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n");
        return false;
    }
    // we set anything masked to NAN except if CONV.POOR is the only bit set
    // First check the old value
    psU32 convPoor = psMetadataLookupU32(&status, masks, "POOR.WARP");
    if (!status) {
        convPoor = psMetadataLookupU32(&status, masks, "CONV.POOR");
        if (!status) {
            psError(PM_ERR_CONFIG, false, "failed to lookup mask value for CONV.POOR in recipes\n");
            return false;
        }
    }
    psU32 maskMask = ~convPoor;

    double exciseValue;
    if (image->type.type == PS_TYPE_U16) {
        exciseValue = 0xffff;
    } else if (image->type.type == PS_TYPE_F32) {
        exciseValue = NAN;
    } else {
         psError(PS_ERR_PROGRAMMING, true, "unexpected image type: %d\n", image->type.type);
        return false;
    }
    long numExcised = 0;
    for (int y=0; y<image->numRows; y++) {
        for (int x=0; x<image->numCols; x++) {
            psU16 maskVal = psImageGet(mask, x, y);
            if (maskVal & maskMask) {
                numExcised++;
                psImageSet(image, x, y, exciseValue);
                if (variance) {
                    psImageSet(variance, x, y, exciseValue);
                }
            }
        }
    }
    fprintf(stderr, "excised %ld masked pixels\n", numExcised);

    return true;
}
