#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>

#define NUM 8
#define DET_SIZE 2000
#define SKY_SIZE 1000
#define SIZE 1024
#define OUTPUT_ROOT "test"
#define SCALE 0.654321
#define ROT M_PI_2
#define INTERPOLATION PS_INTERPOLATE_LANCZOS3
#define OFFSET 16
#define SMOOTH_SIGMA 6.54321
#define SMOOTH_N_SIGMA 2.0
#define DUAL_KERNEL "sub.subkernel"

static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 };

static const char *rootNames[] = { "o5298g0209o.fake.warp", "o5298g0210o.fake.warp", "o5298g0211o.fake.warp",
                                   "o5298g0212o.fake.warp", "o5298g0213o.fake.warp", "o5298g0214o.fake.warp",
                                   "o5298g0215o.fake.warp", "o5298g0216o.fake.warp" };
static const char *kernels[] = { "test.5.kernel", "test.11.kernel", "test.17.kernel", "test.23.kernel",
                                 "test.29.kernel", "test.35.kernel", "test.41.kernel", "test.47.kernel" };


void writeImage(const psImage *image, const char *suffix)
{
    psString name = NULL;
    psStringAppend(&name, "%s.%s", OUTPUT_ROOT, suffix);
    psFits *fits = psFitsOpen(name, "w");
    psFree(name);
    psFitsWriteImage(fits, NULL, image, 0, NULL);
    psFitsClose(fits);
}

#if 0
psImage *readImage(const char *name)
{
    psFits *fits = psFitsOpen(name, "r");
    psImage *image = psFitsReadImage(fits, psRegionSet(0,DET_SIZE,0,DET_SIZE), 0);
    psFitsClose(fits);
    return image;
}

psImage *readMask(const char *name)
{
    psFits *fits = psFitsOpen(name, "r");
    psImage *image = psFitsReadImage(fits, psRegionSet(0,DET_SIZE,0,DET_SIZE), 0);
    if (image->type.type != PS_TYPE_IMAGE_MASK) {
        psImage *temp = psImageCopy(NULL, image, PS_TYPE_IMAGE_MASK);
        psFree(image);
        image = temp;
    }
    psFitsClose(fits);
    return image;
}

psImage *readVariance(psKernel **covar, const char *name)
{
    psFits *fits = psFitsOpen(name, "r");
    psImage *variance = psFitsReadImage(fits, psRegionSet(0,DET_SIZE,0,DET_SIZE), 0);
    if (psFitsMoveExtNum(fits, 1, true)) {
        psImage *covariance = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0);
        *covar = psKernelAllocFromImage(covariance, covariance->numCols / 2, covariance->numRows / 2);
        psFree(covariance);
    } else {
        psErrorClear();
    }
    psFitsClose(fits);
    return variance;
}
#endif

float meanVar(const psImage *var, const psImage *mask, const psKernel *covar)
{
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    psImageBackground(stats, NULL, var, mask, 0xFF, rng);
    float variance = stats->sampleMedian * psImageCovarianceFactor(covar);
    psFree(stats);
    psFree(rng);
    return variance;
}


float signoise(const psImage *image, const psImage *mask, const psImage *var, const psKernel *covar)
{
    float varFactor = psImageCovarianceFactor(covar);
    psImage *sn = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
    for (int y = 0; y < image->numRows; y++) {
        for (int x = 0; x < image->numCols; x++) {
            sn->data.F32[y][x] = image->data.F32[y][x] / sqrtf(var->data.F32[y][x] * varFactor);
        }
    }
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV);
    psImageBackground(stats, NULL, sn, mask, 0xFF, rng);
    float noise = stats->robustStdev;
    psFree(stats);
    psFree(rng);
    psFree(sn);
    return noise;
}

void phot(psImage *image, psImage *mask, psImage *variance, psKernel *covar)
{
    psImage *smoothImage = psImageCopy(NULL, image, PS_TYPE_F32);
    psImageSmoothMask(smoothImage, smoothImage, mask, 0xFF, SMOOTH_SIGMA, SMOOTH_N_SIGMA, 0.1);
    psImage *smoothVariance = psImageCopy(NULL, variance, PS_TYPE_F32);
    psImageSmoothMask(smoothVariance, smoothVariance, mask, 0xFF,
                      SMOOTH_SIGMA * M_SQRT1_2, SMOOTH_N_SIGMA, 0.1);
    int extent = SMOOTH_SIGMA * SMOOTH_N_SIGMA + 0.5;
    psImage *smoothMask = psImageConvolveMask(NULL, mask, 0xFF, 0xFF, -extent, extent, -extent, extent);

    psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
    psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar);
    psFree(kernel);
    psImageCovarianceTransfer(smoothVariance, smoothCovar);

    fprintf(stderr, "Phot: S/N: %f Covar: %f Var: %f\n",
            signoise(smoothImage, smoothMask, smoothVariance, smoothCovar),
            psImageCovarianceFactor(smoothCovar),
            meanVar(smoothVariance, smoothMask, smoothCovar));

    psFree(smoothImage);
    psFree(smoothMask);
    psFree(smoothVariance);
    psFree(smoothCovar);
}

void stack(psImage **stackImage, psImage **stackMask, psImage **stackVariance, psKernel **stackCovar,
           psArray *readouts)
{
    int num = readouts->n;
    psArray *covars = psArrayAlloc(num);
    psVector *weights = psVectorAlloc(num, PS_TYPE_F32);
    for (int i = 0; i < num; i++) {
        pmReadout *ro = readouts->data[i];
        weights->data.F32[i] = 1.0 / meanVar(ro->variance, ro->mask, ro->covariance);
        covars->data[i] = psMemIncrRefCounter(ro->covariance);
    }

    *stackImage = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
    *stackMask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
    *stackVariance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
    for (int y = 0; y < SKY_SIZE; y++) {
        for (int x = 0; x < SKY_SIZE; x++) {
            double sumValueWeight = 0.0;
            double sumWeight = 0.0;
            double sumVarianceWeight = 0.0;
            for (int i = 0; i < readouts->n; i++) {
                pmReadout *ro = readouts->data[i];
                sumValueWeight += ro->image->data.F32[y][x] * weights->data.F32[i];
                sumWeight += weights->data.F32[i];
                sumVarianceWeight += ro->variance->data.F32[y][x] * PS_SQR(weights->data.F32[i]);
            }
            (*stackImage)->data.F32[y][x] = sumValueWeight / sumWeight;
            (*stackVariance)->data.F32[y][x] = sumVarianceWeight / PS_SQR(sumWeight);
            (*stackMask)->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = (isfinite((*stackImage)->data.F32[y][x]) && isfinite((*stackVariance)->data.F32[y][x])) ? 0x00 : 0xFF;
        }
    }

    *stackCovar = psImageCovarianceAverageWeighted(covars, weights);
    psFree(weights);
    psFree(covars);
}

void diff(psImage **diffImage, psImage **diffMask, psImage **diffVariance, psKernel **diffCovar,
          const psImage *image1, const psImage *mask1, const psImage *variance1, const psKernel *covar1,
          const psImage *image2, const psImage *mask2, const psImage *variance2, const psKernel *covar2)
{
    *diffImage = (psImage*)psBinaryOp(NULL, (psPtr)image1, "-", (psPtr)image2);
    *diffMask = (psImage*)psBinaryOp(NULL, (psPtr)mask1, "|", (psPtr)mask2);
    *diffVariance = (psImage*)psBinaryOp(NULL, (psPtr)variance1, "+", (psPtr)variance2);
    psArray *covars = psArrayAlloc(2);
    covars->data[0] = psMemIncrRefCounter((psPtr)covar1);
    covars->data[1] = psMemIncrRefCounter((psPtr)covar2);
    psVector *weights = psVectorAlloc(2, PS_TYPE_F32);
    weights->data.F32[0] = meanVar(variance1, mask1, covar1);
    weights->data.F32[1] = meanVar(variance2, mask2, covar2);
    *diffCovar = psImageCovarianceAverageWeighted(covars, weights);
    psFree(covars);
    psFree(weights);
}


void warp(psImage **outImage, psImage **outMask, psImage **outVariance, psKernel **outCovar,
          const psImage *inImage, const psImage *inMask, const psImage *inVariance, const psKernel *inCovar)
{
    psImageInterpolation *interp = psImageInterpolationAlloc(INTERPOLATION, inImage,
                                                             inVariance, inMask, 0xFF,
                                                             NAN, NAN, 0xF0, 0x0F, 0.1,
                                                             1000);
    *outImage = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
    *outMask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
    *outVariance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);

    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    float xOffset = psRandomUniform(rng) * 2 * OFFSET - OFFSET;
    float yOffset = psRandomUniform(rng) * 2 * OFFSET - OFFSET;
    for (int y = 0; y < SKY_SIZE; y++) {
        for (int x = 0; x < SKY_SIZE; x++) {
            float xIn = SCALE * ((x - SKY_SIZE/2) * cos(ROT) + (y - SKY_SIZE/2) * sin(ROT)) + DET_SIZE / 2 + xOffset;
            float yIn = SCALE * ((x - SKY_SIZE/2) * -sin(ROT) + (y - SKY_SIZE/2) * cos(ROT)) + DET_SIZE / 2 + yOffset;

            double img, var;
            psImageMaskType msk;
            psImageInterpolate(&img, &var, &msk, xIn, yIn, interp);

            (*outImage)->data.F32[y][x] = img;
            (*outVariance)->data.F32[y][x] = var;
            (*outMask)->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = isfinite(img) || isfinite(var) ? msk : 0xFF;
        }
    }

    psArray *covariances = psArrayAlloc(1000);
    for (int i = 0; i < covariances->n; i++) {
        psKernel *kernel = psImageInterpolationKernel(psRandomUniform(rng), psRandomUniform(rng),
                                                      INTERPOLATION);
        covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
        psFree(kernel);
    }
    psFree(rng);
    psKernel *avgCovar = psImageCovarianceAverage(covariances);
    psFree(covariances);
    *outCovar = psImageCovarianceScale(avgCovar, SCALE);
    psFree(avgCovar);
}

int main(int argc, char *argv[])
{
    psArray *readouts = psArrayAlloc(NUM);
    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    for (int i = 0; i < NUM; i++) {
        psImage *inImage = psImageAlloc(DET_SIZE, DET_SIZE, PS_TYPE_F32);
        psImage *inMask = psImageAlloc(DET_SIZE, DET_SIZE, PS_TYPE_IMAGE_MASK);
        psImage *inVariance = psImageAlloc(DET_SIZE, DET_SIZE, PS_TYPE_F32);
        psKernel *inCovar = NULL;

        float noise = sqrtf(variances[i]);
        for (int y = 0; y < DET_SIZE; y++) {
            for (int x = 0; x < DET_SIZE; x++) {
                inImage->data.F32[y][x] = psRandomGaussian(rng) * noise;
                inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
                inVariance->data.F32[y][x] = variances[i];
            }
        }

#if 0
        psString imageName = NULL, maskName = NULL, varName = NULL;
        psStringAppend(&imageName, "%s.fits", rootNames[i]);
        psStringAppend(&maskName, "%s.mask.fits", rootNames[i]);
        psStringAppend(&varName, "%s.wt.fits", rootNames[i]);
        inImage = readImage(imageName);
        inMask = readMask(maskName);
        inVariance = readVariance(&inCovar, varName);
        for (int y = 0; y < DET_SIZE; y++) {
            for (int x = 0; x < DET_SIZE; x++) {
                if (!isfinite(inImage->data.F32[y][x]) || !isfinite(inVariance->data.F32[y][x])) {
                    inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0xFF;
                }
            }
        }
#endif

        psImageCovarianceTransfer(inVariance, inCovar);

        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n",
                i,
                signoise(inImage, inMask, inVariance, inCovar),
                psImageCovarianceFactor(inCovar),
                meanVar(inVariance, inMask, inCovar));

        psImage *warpImage = NULL, *warpMask = NULL, *warpVariance = NULL;
        psKernel *warpCovar = NULL;
        warp(&warpImage, &warpMask, &warpVariance, &warpCovar, inImage, inMask, inVariance, inCovar);
        psFree(inImage);
        psFree(inMask);
        psFree(inVariance);
        psFree(inCovar);

        psImageCovarianceTransfer(warpVariance, warpCovar);

        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n",
                i,
                signoise(warpImage, warpMask, warpVariance, warpCovar),
                psImageCovarianceFactor(warpCovar),
                meanVar(warpVariance, warpMask, warpCovar));

        phot(warpImage, warpMask, warpVariance, warpCovar);

        pmReadout *ro = pmReadoutAlloc(NULL);
        ro->image = warpImage;
        ro->mask = warpMask;
        ro->variance = warpVariance;
        ro->covariance = warpCovar;

        psFits *fits = psFitsOpen(kernels[i], "r");
        pmReadoutReadSubtractionKernels(ro, fits);
        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
        kernels->xMax = SKY_SIZE;
        kernels->yMax = SKY_SIZE;
        psFitsClose(fits);

        pmReadout *conv = pmReadoutAlloc(NULL);
        conv->image = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        conv->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
        conv->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.001,
                                       0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
            psErrorStackPrint(stderr, "Error:");
            exit(1);
        }
        psFree(ro);

        psImageCovarianceTransfer(conv->variance, conv->covariance);

        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n",
                i,
                signoise(conv->image, conv->mask, conv->variance, conv->covariance),
                psImageCovarianceFactor(conv->covariance),
                meanVar(conv->variance, conv->mask, conv->covariance));

        phot(conv->image, conv->mask, conv->variance, conv->covariance);
        readouts->data[i] = conv;
    }

    {
        pmReadout *minuend = readouts->data[0];
        pmReadout *subtrahend = readouts->data[readouts->n - 1];
        psImage *diffImage = NULL, *diffMask = NULL, *diffVariance = NULL;
        psKernel *diffCovar = NULL;
        diff(&diffImage, &diffMask, &diffVariance, &diffCovar,
             minuend->image, minuend->mask, minuend->variance, minuend->covariance,
             subtrahend->image, subtrahend->mask, subtrahend->variance, subtrahend->covariance);
        psImageCovarianceTransfer(diffVariance, diffCovar);
        fprintf(stderr, "WW Diff Image: S/N: %f Covar: %f Var: %f\n",
                signoise(diffImage, diffMask, diffVariance, diffCovar),
                psImageCovarianceFactor(diffCovar),
                meanVar(diffVariance, diffMask, diffCovar));

        writeImage(diffImage, "wwdiff.image.fits");
        writeImage(diffMask, "wwdiff.mask.fits");
        writeImage(diffVariance, "wwdiff.variance.fits");
        writeImage(diffCovar->image, "wwdiff.covar.fits");

        psFree(diffImage);
        psFree(diffMask);
        psFree(diffVariance);
        psFree(diffCovar);
    }

    // stack-stack diff
    {
        psArray *readouts1 = psArrayAllocEmpty(NUM);
        psArray *readouts2 = psArrayAllocEmpty(NUM);
        for (int i = 0; i < NUM; i++) {
            if (i <= NUM / 2) {
                psArrayAdd(readouts1, 0, readouts->data[i]);
            } else {
                psArrayAdd(readouts2, 0, readouts->data[i]);
            }
        }

        psImage *stackImage1 = NULL, *stackMask1 = NULL, *stackVariance1 = NULL;
        psKernel *stackCovar1 = NULL;
        psImage *stackImage2 = NULL, *stackMask2 = NULL, *stackVariance2 = NULL;
        psKernel *stackCovar2 = NULL;
        stack(&stackImage1, &stackMask1, &stackVariance1, &stackCovar1, readouts1);
        stack(&stackImage2, &stackMask2, &stackVariance2, &stackCovar2, readouts2);
        psFree(readouts1);
        psFree(readouts2);
        psImageCovarianceTransfer(stackVariance1, stackCovar1);
        psImageCovarianceTransfer(stackVariance2, stackCovar2);
        fprintf(stderr, "Stack 1: S/N: %f Covar: %f Var: %f\n",
                signoise(stackImage1, stackMask1, stackVariance1, stackCovar1),
                psImageCovarianceFactor(stackCovar1),
                meanVar(stackVariance1, stackMask1, stackCovar1));
        fprintf(stderr, "Stack 2: S/N: %f Covar: %f Var: %f\n",
                signoise(stackImage2, stackMask2, stackVariance2, stackCovar2),
                psImageCovarianceFactor(stackCovar2),
                meanVar(stackVariance2, stackMask2, stackCovar2));
        phot(stackImage1, stackMask1, stackVariance1, stackCovar1);
        phot(stackImage2, stackMask2, stackVariance2, stackCovar2);

        pmReadout *stack1 = pmReadoutAlloc(NULL);
        pmReadout *stack2 = pmReadoutAlloc(NULL);
        stack1->image = stackImage1;
        stack1->mask = stackMask1;
        stack1->variance = stackVariance1;
        stack1->covariance = stackCovar1;
        stack2->image = stackImage2;
        stack2->mask = stackMask2;
        stack2->variance = stackVariance2;
        stack2->covariance = stackCovar2;

        psFits *fits = psFitsOpen(DUAL_KERNEL, "r");
        pmReadoutReadSubtractionKernels(stack1, fits);
        psFitsClose(fits);
        pmSubtractionKernels *kernel = psMetadataLookupPtr(NULL, stack1->analysis,
                                                           PM_SUBTRACTION_ANALYSIS_KERNEL);
        kernel->xMax = SKY_SIZE;
        kernel->yMax = SKY_SIZE;

        pmReadout *convStack1 = pmReadoutAlloc(NULL);
        pmReadout *convStack2 = pmReadoutAlloc(NULL);
        convStack1->image = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        convStack1->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
        convStack1->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        convStack2->image = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        convStack2->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
        convStack2->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);

        if (!pmSubtractionMatchPrecalc(convStack1, convStack2, stack1, stack2, stack1->analysis,
                                       32, 0.0, 0.001, 0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
            psErrorStackPrint(stderr, "Error:");
            exit(1);
        }
        psFree(stack1);
        psFree(stack2);

        psImageCovarianceTransfer(convStack1->variance, convStack1->covariance);
        psImageCovarianceTransfer(convStack2->variance, convStack2->covariance);

        fprintf(stderr, "Conv stack 1: S/N: %f Covar: %f Var: %f\n",
                signoise(convStack1->image, convStack1->mask, convStack1->variance, convStack1->covariance),
                psImageCovarianceFactor(convStack1->covariance),
                meanVar(convStack1->variance, convStack1->mask, convStack1->covariance));
        fprintf(stderr, "Conv stack 2: S/N: %f Covar: %f Var: %f\n",
                signoise(convStack2->image, convStack2->mask, convStack2->variance, convStack2->covariance),
                psImageCovarianceFactor(convStack2->covariance),
                meanVar(convStack2->variance, convStack2->mask, convStack2->covariance));

        psImage *diffImage = NULL, *diffMask = NULL, *diffVariance = NULL;
        psKernel *diffCovar = NULL;
        diff(&diffImage, &diffMask, &diffVariance, &diffCovar,
             convStack1->image, convStack1->mask, convStack1->variance, convStack1->covariance,
             convStack2->image, convStack2->mask, convStack2->variance, convStack2->covariance);
        psImageCovarianceTransfer(diffVariance, diffCovar);
        fprintf(stderr, "SS Diff Image: S/N: %f Covar: %f Var: %f\n",
                signoise(diffImage, diffMask, diffVariance, diffCovar),
                psImageCovarianceFactor(diffCovar),
                meanVar(diffVariance, diffMask, diffCovar));

        writeImage(diffImage, "ssdiff.image.fits");
        writeImage(diffMask, "ssdiff.mask.fits");
        writeImage(diffVariance, "ssdiff.variance.fits");
        writeImage(diffCovar->image, "ssdiff.covar.fits");

        phot(diffImage, diffMask, diffVariance, diffCovar);

        psFree(diffImage);
        psFree(diffMask);
        psFree(diffVariance);
        psFree(diffCovar);

        psFree(convStack1);
        psFree(convStack2);
    }

    // Full stack and warp-stack diff
    {
        psImage *stackImage = NULL, *stackMask = NULL, *stackVariance = NULL;
        psKernel *stackCovar = NULL;
        stack(&stackImage, &stackMask, &stackVariance, &stackCovar, readouts);
        psImageCovarianceTransfer(stackVariance, stackCovar);

        writeImage(stackImage, "stack.image.fits");
        writeImage(stackMask, "stack.mask.fits");
        writeImage(stackVariance, "stack.variance.fits");
        writeImage(stackCovar->image, "stack.covar.fits");

        fprintf(stderr, "Stack: S/N: %f Covar: %f Var: %f\n",
                signoise(stackImage, stackMask, stackVariance, stackCovar),
                psImageCovarianceFactor(stackCovar),
                meanVar(stackVariance, stackMask, stackCovar));

        phot(stackImage, stackMask, stackVariance, stackCovar);

        pmReadout *stack = pmReadoutAlloc(NULL);
        stack->image = stackImage;
        stack->mask = stackMask;
        stack->variance = stackVariance;
        stack->covariance = stackCovar;

        psFits *fits = psFitsOpen(DUAL_KERNEL, "r");
        pmReadoutReadSubtractionKernels(stack, fits);
        psFitsClose(fits);
        pmSubtractionKernels *kernel = psMetadataLookupPtr(NULL, stack->analysis,
                                                           PM_SUBTRACTION_ANALYSIS_KERNEL);
        kernel->xMax = SKY_SIZE;
        kernel->yMax = SKY_SIZE;

        pmReadout *warp = readouts->data[readouts->n - 1];
        pmReadout *convStack = pmReadoutAlloc(NULL);
        pmReadout *convWarp = pmReadoutAlloc(NULL);
        convStack->image = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        convStack->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
        convStack->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        convWarp->image = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
        convWarp->mask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
        convWarp->variance = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);

        if (!pmSubtractionMatchPrecalc(convStack, convWarp, stack, warp, stack->analysis,
                                       32, 0.0, 0.001, 0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
            psErrorStackPrint(stderr, "Error:");
            exit(1);
        }
        psFree(stack);
        psFree(readouts);

        psImageCovarianceTransfer(convStack->variance, convStack->covariance);
        psImageCovarianceTransfer(convWarp->variance, convWarp->covariance);

        fprintf(stderr, "Conv stack: S/N: %f Covar: %f Var: %f\n",
                signoise(convStack->image, convStack->mask, convStack->variance, convStack->covariance),
                psImageCovarianceFactor(convStack->covariance),
                meanVar(convStack->variance, convStack->mask, convStack->covariance));
        fprintf(stderr, "Conv warp: S/N: %f Covar: %f Var: %f\n",
                signoise(convWarp->image, convWarp->mask, convWarp->variance, convWarp->covariance),
                psImageCovarianceFactor(convWarp->covariance),
                meanVar(convWarp->variance, convWarp->mask, convWarp->covariance));

        psImage *diffImage = NULL, *diffMask = NULL, *diffVariance = NULL;
        psKernel *diffCovar = NULL;
        diff(&diffImage, &diffMask, &diffVariance, &diffCovar,
             convStack->image, convStack->mask, convStack->variance, convStack->covariance,
             convWarp->image, convWarp->mask, convWarp->variance, convWarp->covariance);
        psFree(convStack);
        psFree(convWarp);
        psImageCovarianceTransfer(diffVariance, diffCovar);
        fprintf(stderr, "WS Diff Image: S/N: %f Covar: %f Var: %f\n",
                signoise(diffImage, diffMask, diffVariance, diffCovar),
                psImageCovarianceFactor(diffCovar),
                meanVar(diffVariance, diffMask, diffCovar));

        writeImage(diffImage, "wsdiff.image.fits");
        writeImage(diffMask, "wsdiff.mask.fits");
        writeImage(diffVariance, "wsdiff.variance.fits");
        writeImage(diffCovar->image, "wsdiff.covar.fits");

        phot(diffImage, diffMask, diffVariance, diffCovar);

        psFree(diffImage);
        psFree(diffMask);
        psFree(diffVariance);
        psFree(diffCovar);
    }

}
