Changeset 27791

Timestamp:

Apr 28, 2010, 12:13:30 PM (16 years ago)

Author:

Paul Price

Message:

Adding more simulations. The entire pipeline is more or less simulated, measuring the noise model at each stage.

Location:

branches/pap/archive/variance_covariance

Files:

• simulate.c (moved) (moved from branches/pap/archive/variance_covariance/stack.c ) (5 diffs)
• sub.subkernel (added)

branches/pap/archive/variance_covariance/simulate.c

@@ -7 +7 @@
 #define SKY_SIZE 1000
 #define SIZE 1024
-#define OUTPUT_ROOT "stack"
+#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 };
@@ -21 +24 @@
                                  "test.29.kernel", "test.35.kernel", "test.41.kernel", "test.47.kernel" };
 
-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);
-}
 
 void writeImage(const psImage *image, const char *suffix)
@@ -143 +104 @@
 }
 
+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);
-    psArray *covars = psArrayAlloc(NUM);
-    psVector *weights = psVectorAlloc(NUM, PS_TYPE_F32);
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     for (int i = 0; i < NUM; i++) {
@@ -206 +289 @@
                 meanVar(warpVariance, warpMask, warpCovar));
 
+        phot(warpImage, warpMask, warpVariance, warpCovar);
+
         pmReadout *ro = pmReadoutAlloc(NULL);
         ro->image = warpImage;
@@ -238 +323 @@
                 meanVar(conv->variance, conv->mask, conv->covariance));
 
-        weights->data.F32[i] = 1.0 / meanVar(conv->variance, conv->mask, conv->covariance);
-
-        covars->data[i] = psMemIncrRefCounter(conv->covariance);
+        phot(conv->image, conv->mask, conv->variance, conv->covariance);
         readouts->data[i] = conv;
     }
 
-    psImage *stackImage = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_F32);
-    psImage *stackMask = psImageAlloc(SKY_SIZE, SKY_SIZE, PS_TYPE_IMAGE_MASK);
-    psImage *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 < NUM; 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]);
+    {
+        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]);
             }
-            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;
-        }
-    }
-
-    psKernel *stackCovar = psImageCovarianceAverageWeighted(covars, weights);
-    psFree(weights);
-    psFree(covars);
-    psFree(readouts);
-
-    psImageCovarianceTransfer(stackVariance, stackCovar);
-
-    writeImage(stackImage, "image.fits");
-    writeImage(stackMask, "mask.fits");
-    writeImage(stackVariance, "variance.fits");
-    writeImage(stackCovar->image, "covar.fits");
-
-    fprintf(stderr, "Stack: S/N: %f Covar: %f Var: %f\n",
-            signoise(stackImage, stackMask, stackVariance, stackCovar),
-            psImageCovarianceFactor(stackCovar),
-            meanVar(stackVariance, stackMask, stackCovar));
-
-    float smoothSigma = 6.680111;
-    float smoothNsigma = 2.000000;
-
-    psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
-
-    psKernel *smoothCovar = psImageCovarianceCalculate(kernel, stackCovar);
-    psImage *smoothImage = psImageCopy(NULL, stackImage, PS_TYPE_F32);
-    psImageSmoothMask(smoothImage, smoothImage, stackMask, 0xFF,
-                      smoothSigma, smoothNsigma, 0.1);
-    psImage *smoothVariance = psImageCopy(NULL, stackVariance, PS_TYPE_F32);
-    psImageSmoothMask(smoothVariance, smoothVariance, stackMask, 0xFF,
-                      smoothSigma * M_SQRT1_2, smoothNsigma, 0.1);
-    int extent = smoothSigma * smoothNsigma + 0.5;
-    psImage *smoothMask = psImageConvolveMask(NULL, stackMask, 0xFF, 0xFF, -extent, extent, -extent, extent);
-
-    psImageCovarianceTransfer(smoothVariance, smoothCovar);
-
-    fprintf(stderr, "Smoothed: S/N: %f Covar: %f Var: %f\n",
-            signoise(smoothImage, smoothMask, smoothVariance, smoothCovar),
-            psImageCovarianceFactor(smoothCovar),
-            meanVar(smoothVariance, smoothMask, smoothCovar));
-
-    psFree(stackImage);
-    psFree(stackMask);
-    psFree(stackVariance);
-    psFree(stackCovar);
-
-    psFree(smoothImage);
-    psFree(smoothMask);
-    psFree(smoothVariance);
-    psFree(smoothCovar);
-}
+        }
+
+        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);
+    }
+
+}