Changeset 28159

Timestamp:

May 28, 2010, 12:04:26 PM (16 years ago)

Author:

Paul Price

Message:

Bug fixes, getting ready for publication. Added histogram S/N ouput.

File:

• trunk/archive/noise_model/simulate.c (modified) (10 diffs)

trunk/archive/noise_model/simulate.c

@@ -10 +10 @@
 #define SCALE 0.654321
 #define ROT M_PI_2
-#define INTERPOLATION PS_INTERPOLATE_LANCZOS3
+#define INTERPOLATION PS_INTERPOLATE_LANCZOS4
 #define OFFSET 16
 #define SMOOTH_SIGMA 6.54321
 #define SMOOTH_N_SIGMA 2.0
 #define DUAL_KERNEL "sub.subkernel"
+#define WARP_NUM 10000
 
 static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 };
@@ -94 +95 @@
         }
     }
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV);
-    psImageBackground(stats, NULL, sn, mask, 0xFF, rng);
-    float noise = stats->robustStdev;
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_STDEV);
+    psImageStats(stats, sn, mask, 0xFF);
+    float noise = stats->sampleStdev;
     psFree(stats);
-    psFree(rng);
+
+    psHistogram *hist = psHistogramAlloc(-5, +5, 1001);
+    psVector *data = psVectorAlloc(sn->numCols * sn->numRows, PS_TYPE_F32);
+    psVector *dataMask = psVectorAlloc(sn->numCols * sn->numRows, PS_TYPE_VECTOR_MASK);
+    psVectorInit(dataMask, 0);
+    long num = 0;
+    for (int y = 0, i = 0; y < sn->numRows; y++) {
+        for (int x = 0; x < sn->numCols; x++, i++) {
+            data->data.F32[i] = sn->data.F32[y][x];
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) {
+                dataMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF;
+            } else {
+                num++;
+            }
+        }
+    }
+    psVectorHistogram(hist, data, NULL, dataMask, 0xFF);
+    psFree(data);
+    psFree(dataMask);
     psFree(sn);
+
+    static int number = 0;
+    psString name = NULL;
+    fprintf(stderr, "Writing histogram %d\n", number);
+    psStringAppend(&name, "hist_%d.dat", number++);
+    FILE *file = fopen(name, "w");
+    psFree(name);
+    fprintf(stderr, "# Sig Frac\n");
+    for (int i = 0; i < hist->bounds->n - 1; i++) {
+        fprintf(file, "%f %f\n",
+                0.5 * (hist->bounds->data.F32[i] + hist->bounds->data.F32[i+1]),
+                hist->nums->data.F32[i] / num);
+    }
+    fclose(file);
+    psFree(hist);
+
     return noise;
 }
 
+void convolve(psImage **smoothImage, psImage **smoothMask, psImage **smoothVariance, psKernel **smoothCovar,
+              const psImage *image, const psImage *mask, const psImage *variance, const psKernel *covar,
+              const psKernel *kernel)
+{
+    psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
+    double sum = 0.0, sum2 = 0.0;
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sum += kernel->kernel[y][x];
+            sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            kernel2->kernel[y][x] /= sum2;
+        }
+    }
+    *smoothImage = psImageConvolveDirect(NULL, image, kernel);
+    *smoothVariance = psImageConvolveDirect(NULL, variance, kernel2);
+    *smoothCovar = psImageCovarianceCalculate(kernel, covar);
+    *smoothMask = psImageCopy(NULL, mask, PS_TYPE_IMAGE_MASK); // Cheating
+    psFree(kernel2);
+}
+
+
 void phot(psImage *image, psImage *mask, psImage *variance, psKernel *covar)
 {
+#if 1
     psImage *smoothImage = psImageCopy(NULL, image, PS_TYPE_F32);
+    psImage *smoothVariance = psImageCopy(NULL, variance, 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);
+                      SMOOTH_SIGMA * M_SQRT1_2, SMOOTH_N_SIGMA / M_SQRT1_2, 0.1);
     int extent = SMOOTH_SIGMA * SMOOTH_N_SIGMA + 0.5;
     psImage *smoothMask = psImageConvolveMask(NULL, mask, 0xFF, 0xFF, -extent, extent, -extent, extent);
@@ -118 +178 @@
     psFree(kernel);
     psImageCovarianceTransfer(smoothVariance, smoothCovar);
+#else
+    psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
+    psKernel *kernel2 = psKernelAlloc(kernel->xMin, kernel->xMax, kernel->yMin, kernel->yMax);
+    double sum = 0.0, sum2 = 0.0;
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sum += kernel->kernel[y][x];
+            sum2 += kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    fprintf(stderr, "Kernel sum: %f %f\n", sum, sum2);
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            kernel2->kernel[y][x] /= sum2;
+        }
+    }
+    psImage *smoothImage = psImageConvolveDirect(NULL, image, kernel);
+    psImage *smoothVariance = psImageConvolveDirect(NULL, variance, kernel2);
+    psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar);
+    psImage *smoothMask = psImageCopy(NULL, mask, PS_TYPE_IMAGE_MASK);
+#endif
 
     fprintf(stderr, "Phot: S/N: %f Covar: %f Var: %f\n",
@@ -200 +281 @@
     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;
+    for (int y = 0, index = 0; y < SKY_SIZE; y++) {
+        float dy = y - SKY_SIZE/2;
+        for (int x = 0; x < SKY_SIZE; x++, index++) {
+            float dx = x - SKY_SIZE/2;
+            float xIn = SCALE * (dx * cos(ROT) + dy * sin(ROT)) + DET_SIZE / 2 + xOffset;
+            float yIn = SCALE * (dx * -sin(ROT) + dy * cos(ROT)) + DET_SIZE / 2 + yOffset;
+
+#if 0
+            xIn += 0.12345e-6 * PS_SQR(dx);
+            yIn += 0.12345e-6 * PS_SQR(dy);
+#endif
+
+            double img = NAN, var = NAN;
+            psImageMaskType msk = 0;
             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++) {
+            (*outMask)->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = isfinite(img) && isfinite(var) ? msk : 0xFF;
+        }
+    }
+
+    psArray *covariances = psArrayAlloc(WARP_NUM);
+    psVector *factors = psVectorAlloc(WARP_NUM, PS_TYPE_F32);
+    double mean = 0.0;
+    for (int i = 0; i < WARP_NUM; i++) {
         psKernel *kernel = psImageInterpolationKernel(psRandomUniform(rng), psRandomUniform(rng),
                                                       INTERPOLATION);
-        covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
+        psKernel *covar = covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
         psFree(kernel);
+        mean += factors->data.F32[i] = psImageCovarianceFactor(covar);
     }
     psFree(rng);
     psKernel *avgCovar = psImageCovarianceAverage(covariances);
     psFree(covariances);
+
     *outCovar = psImageCovarianceScale(avgCovar, SCALE);
     psFree(avgCovar);
+
+    mean /= WARP_NUM;
+
+    double stdev = 0.0;
+    for (int i = 0; i < WARP_NUM; i++) {
+        stdev += PS_SQR(factors->data.F32[i] - mean);
+    }
+    stdev = sqrt(stdev/(WARP_NUM-1));
+    fprintf(stderr, "Warp covariance mean: %f stdev: %f\n", mean, stdev);
+    psFree(factors);
 }
 
@@ -265 +367 @@
 #endif
 
-        psImageCovarianceTransfer(inVariance, inCovar);
+        //        psImageCovarianceTransfer(inVariance, inCovar);
+
+        {
+            psString imageName = NULL, maskName = NULL, varName = NULL;
+            psStringAppend(&imageName, "input.image.%d.fits", i);
+            psStringAppend(&maskName, "input.mask.%d.fits", i);
+            psStringAppend(&varName, "input.var.%d.fits", i);
+            writeImage(inImage, imageName);
+            writeImage(inMask, maskName);
+            writeImage(inVariance, varName);
+            psFree(imageName);
+            psFree(maskName);
+            psFree(varName);
+        }
 
         fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n",
@@ -272 +387 @@
                 psImageCovarianceFactor(inCovar),
                 meanVar(inVariance, inMask, inCovar));
+
+        phot(inImage, inMask, inVariance, inCovar);
 
         psImage *warpImage = NULL, *warpMask = NULL, *warpVariance = NULL;
@@ -281 +398 @@
         psFree(inCovar);
 
-        psImageCovarianceTransfer(warpVariance, warpCovar);
+        //        psImageCovarianceTransfer(warpVariance, warpCovar);
+
+        {
+            psString imageName = NULL, maskName = NULL, varName = NULL, covarName = NULL;
+            psStringAppend(&imageName, "warp.image.%d.fits", i);
+            psStringAppend(&maskName, "warp.mask.%d.fits", i);
+            psStringAppend(&varName, "warp.var.%d.fits", i);
+            psStringAppend(&covarName, "warp.covar.%d.fits", i);
+            writeImage(warpImage, imageName);
+            writeImage(warpMask, maskName);
+            writeImage(warpVariance, varName);
+            writeImage(warpCovar->image, covarName);
+            psFree(imageName);
+            psFree(maskName);
+            psFree(varName);
+            psFree(covarName);
+        }
 
         fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n",
@@ -300 +433 @@
         pmReadoutReadSubtractionKernels(ro, fits);
         pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
+#if 0
+        for (int i = 0; i < kernels->num; i++) {
+            for (int y = 0, index = i; y < kernels->spatialOrder; y++) {
+                for (int x = 0; x < kernels->spatialOrder - y; x++, index++) {
+                    if (x != 0 && y != 0) {
+                        if (kernels->solution1) {
+                            kernels->solution1->data.F64[index] = 0.0;
+                        }
+                        if (kernels->solution2) {
+                            kernels->solution2->data.F64[index] = 0.0;
+                        }
+                    }
+                }
+            }
+        }
+#endif
+
         kernels->xMax = SKY_SIZE;
         kernels->yMax = SKY_SIZE;
@@ -305 +455 @@
 
         pmReadout *conv = pmReadoutAlloc(NULL);
+#if 1
         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,
+        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 32, 0.0, 0.01,
                                        0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
             psErrorStackPrint(stderr, "Error:");
             exit(1);
         }
+#else
+        {
+            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false);
+            convolve(&conv->image, &conv->mask, &conv->variance, &conv->covariance,
+                     ro->image, ro->mask, ro->variance, ro->covariance,
+                     kernel);
+            psFree(kernel);
+        }
+#endif
         psFree(ro);
 
-        psImageCovarianceTransfer(conv->variance, conv->covariance);
+        //        psImageCovarianceTransfer(conv->variance, conv->covariance);
+
+        {
+            psString imageName = NULL, maskName = NULL, varName = NULL, covarName = NULL;
+            psStringAppend(&imageName, "conv.image.%d.fits", i);
+            psStringAppend(&maskName, "conv.mask.%d.fits", i);
+            psStringAppend(&varName, "conv.var.%d.fits", i);
+            psStringAppend(&covarName, "conv.covar.%d.fits", i);
+            writeImage(conv->image, imageName);
+            writeImage(conv->mask, maskName);
+            writeImage(conv->variance, varName);
+            writeImage(conv->covariance->image, covarName);
+            psFree(imageName);
+            psFree(maskName);
+            psFree(varName);
+            psFree(covarName);
+        }
 
         fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n",
@@ -325 +501 @@
         phot(conv->image, conv->mask, conv->variance, conv->covariance);
         readouts->data[i] = conv;
+
     }