Changeset 28794 for branches/eam_branches/ipp-20100621/archive/noise_model/simulate.c

Timestamp:

Jul 30, 2010, 9:31:50 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

File:

• branches/eam_branches/ipp-20100621/archive/noise_model/simulate.c (modified) (23 diffs)

branches/eam_branches/ipp-20100621/archive/noise_model/simulate.c

@@ -6 +6 @@
 #define DET_SIZE 2000
 #define SKY_SIZE 1000
-#define SIZE 1024
 #define OUTPUT_ROOT "test"
-#define SCALE 0.654321
+#define SCALE 0.7654321
 #define ROT M_PI_2
-#define INTERPOLATION PS_INTERPOLATE_LANCZOS4
+#define INTERPOLATION PS_INTERPOLATE_LANCZOS3
 #define OFFSET 16
 #define SMOOTH_SIGMA 6.54321
@@ -16 +15 @@
 #define DUAL_KERNEL "sub.subkernel"
 #define WARP_NUM 10000
+#define CONV_NUM 1000
 
 static const float variances[] = { 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0 };
@@ -21 +21 @@
 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" };
+#if 1
 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" };
-
+#else
+static const char *kernels[] = { "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel",
+                                 "test.11.kernel", "test.11.kernel", "test.11.kernel", "test.11.kernel" };
+#endif
 
 void writeImage(const psImage *image, const char *suffix)
@@ -76 +81 @@
 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);
+    psImageStats(stats, var, mask, 0xFF);
     float variance = stats->sampleMedian * psImageCovarianceFactor(covar);
     psFree(stats);
-    psFree(rng);
     return variance;
 }
@@ -93 +96 @@
         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);
+            if (!isfinite(sn->data.F32[y][x])) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0xFF;
+            }
         }
     }
@@ -100 +106 @@
     float offset = stats->sampleMean;
     psFree(stats);
+
+    static int number = 0;
+
+    psString snName = NULL;
+    psStringAppend(&snName, "sn.%d.fits", number);
+    writeImage(sn, snName);
+    fprintf(stderr, "Writing S/N image %d\n", number);
 
     psHistogram *hist = psHistogramAlloc(-5, +5, 101);
@@ -121 +134 @@
     psFree(sn);
 
-    static int number = 0;
     psString name = NULL;
     fprintf(stderr, "Writing histogram %d\n", number);
-    psStringAppend(&name, "hist_%d.dat", number++);
+    psStringAppend(&name, "hist_%d.dat", number);
     FILE *file = fopen(name, "w");
     psFree(name);
@@ -136 +148 @@
     psFree(hist);
 
+    number++;
+
     return noise;
 }
@@ -144 +158 @@
 {
     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;
+            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -176 +183 @@
 
     psKernel *kernel = psImageSmoothKernel(SMOOTH_SIGMA, SMOOTH_N_SIGMA); // Kernel used for smoothing
+    double sum2 = 0.0;                                               // Sum of kernel squared
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sum2 += PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    float factor = 1.0 / sum2;
     psKernel *smoothCovar = psImageCovarianceCalculate(kernel, covar);
     psFree(kernel);
-    psImageCovarianceTransfer(smoothVariance, smoothCovar);
+
+    // Apply square root of significance image scaling factor to image
+    psBinaryOp(smoothImage, smoothImage, "*", psScalarAlloc(sqrtf(factor), PS_TYPE_F32));
+
 #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;
+            kernel2->kernel[y][x] = PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -332 +341 @@
 }
 
+float covarianceSum(const psKernel *covar)
+{
+    if (!covar) {
+        return 1.0;
+    }
+    double sum = 0.0;
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            sum += covar->kernel[y][x];
+        }
+    }
+    return sum;
+}
+
+
 int main(int argc, char *argv[])
 {
@@ -383 +407 @@
         }
 
-        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Input image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(inImage, inMask, inVariance, inCovar),
                 psImageCovarianceFactor(inCovar),
-                meanVar(inVariance, inMask, inCovar));
+                meanVar(inVariance, inMask, inCovar),
+                covarianceSum(inCovar));
 
         phot(inImage, inMask, inVariance, inCovar);
@@ -417 +442 @@
         }
 
-        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Warp image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(warpImage, warpMask, warpVariance, warpCovar),
                 psImageCovarianceFactor(warpCovar),
-                meanVar(warpVariance, warpMask, warpCovar));
+                meanVar(warpVariance, warpMask, warpCovar),
+                covarianceSum(warpCovar));
 
         phot(warpImage, warpMask, warpVariance, warpCovar);
@@ -434 +460 @@
         pmReadoutReadSubtractionKernels(ro, fits);
         pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels);
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels);
+#if 1
+//        kernels->solution1->data.F64[normIndex] += 1.0;
+        kernels->solution1->data.F64[bgIndex] = 0.0;
+#endif
+        fprintf(stderr, "Norm: %f BG: %f\n", kernels->solution1->data.F64[normIndex], kernels->solution1->data.F64[bgIndex]);
 #if 0
         for (int i = 0; i < kernels->num; i++) {
@@ -455 +488 @@
         psFitsClose(fits);
 
+#if 0
+        {
+            psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+            psArray *covariances = psArrayAlloc(CONV_NUM);
+            psVector *factors = psVectorAlloc(CONV_NUM, PS_TYPE_F32);
+            double mean = 0.0;
+            for (int i = 0; i < CONV_NUM; i++) {
+                float x, y;
+                p_pmSubtractionPolynomialNormCoords(&x, &y, psRandomUniform(rng) * SKY_SIZE,
+                                                    psRandomUniform(rng) * SKY_SIZE,
+                                                    kernels->xMin, kernels->xMax,
+                                                    kernels->yMin, kernels->yMax);
+                psKernel *kernel = pmSubtractionKernel(kernels, x, y, false);
+                psKernel *covar = covariances->data[i] = psImageCovarianceCalculate(kernel, inCovar);
+                psFree(kernel);
+                mean += factors->data.F32[i] = psImageCovarianceFactor(covar);
+            }
+            psFree(rng);
+            psFree(covariances);
+
+            mean /= WARP_NUM;
+
+            double stdev = 0.0;
+            for (int i = 0; i < CONV_NUM; i++) {
+                stdev += PS_SQR(factors->data.F32[i] - mean);
+            }
+            stdev = sqrt(stdev/(CONV_NUM-1));
+            fprintf(stderr, "Conv covariance mean: %f stdev: %f\n", mean, stdev);
+            psFree(factors);
+        }
+#endif
+
         pmReadout *conv = pmReadoutAlloc(NULL);
 #if 1
@@ -460 +525 @@
         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.01,
+        if (!pmSubtractionMatchPrecalc(NULL, conv, NULL, ro, ro->analysis, 2 * kernels->size + 1, 0.0, 0.01,
                                        0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
             psErrorStackPrint(stderr, "Error:");
@@ -494 +559 @@
         }
 
-        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f\n",
+        fprintf(stderr, "Conv Image %d: S/N: %f Covar: %f Var: %f CovarSum: %f\n",
                 i,
                 signoise(conv->image, conv->mask, conv->variance, conv->covariance),
                 psImageCovarianceFactor(conv->covariance),
-                meanVar(conv->variance, conv->mask, conv->covariance));
+                meanVar(conv->variance, conv->mask, conv->covariance),
+                covarianceSum(conv->covariance));
 
         phot(conv->image, conv->mask, conv->variance, conv->covariance);
         readouts->data[i] = conv;
 
+        //        exit(1);
     }
 
@@ -562 +629 @@
                 psImageCovarianceFactor(stackCovar2),
                 meanVar(stackVariance2, stackMask2, stackCovar2));
+
+        writeImage(stackImage1, "stack1.image.fits");
+        writeImage(stackMask1, "stack1.mask.fits");
+        writeImage(stackVariance1, "stack1.variance.fits");
+        writeImage(stackCovar1->image, "stack1.covar.fits");
+        writeImage(stackImage2, "stack2.image.fits");
+        writeImage(stackMask2, "stack2.mask.fits");
+        writeImage(stackVariance2, "stack2.variance.fits");
+        writeImage(stackCovar2->image, "stack2.covar.fits");
+
         phot(stackImage1, stackMask1, stackVariance1, stackCovar1);
         phot(stackImage2, stackMask2, stackVariance2, stackCovar2);
@@ -584 +661 @@
         kernel->yMax = SKY_SIZE;
 
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernel);
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernel);
+#if 1
+//        kernel->solution1->data.F64[normIndex] += 1.0;
+        kernel->solution1->data.F64[bgIndex] = 0.0;
+#endif
+        fprintf(stderr, "Norm: %f BG: %f\n", kernel->solution1->data.F64[normIndex], kernel->solution1->data.F64[bgIndex]);
+
         pmReadout *convStack1 = pmReadoutAlloc(NULL);
         pmReadout *convStack2 = pmReadoutAlloc(NULL);
@@ -594 +679 @@
 
         if (!pmSubtractionMatchPrecalc(convStack1, convStack2, stack1, stack2, stack1->analysis,
-                                       32, 0.0, 0.001, 0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
+                                       2 * kernel->size + 1, 0.0, 0.001, 0xFF, 0xF0, 0x0F, 0.1, 1.0)) {
             psErrorStackPrint(stderr, "Error:");
             exit(1);
@@ -600 +685 @@
         psFree(stack1);
         psFree(stack2);
+
+        writeImage(convStack1->image, "stack1.conv.image.fits");
+        writeImage(convStack1->mask, "stack1.conv.mask.fits");
+        writeImage(convStack1->variance, "stack1.conv.variance.fits");
+        writeImage(convStack1->covariance->image, "stack1.conv.covar.fits");
+        writeImage(convStack2->image, "stack2.conv.image.fits");
+        writeImage(convStack2->mask, "stack2.conv.mask.fits");
+        writeImage(convStack2->variance, "stack2.conv.variance.fits");
+        writeImage(convStack2->covariance->image, "stack2.conv.covar.fits");
 
         psImageCovarianceTransfer(convStack1->variance, convStack1->covariance);
@@ -624 +718 @@
                 meanVar(diffVariance, diffMask, diffCovar));
 
-        phot(diffImage, diffMask, diffVariance, diffCovar);
-
         writeImage(diffImage, "ssdiff.image.fits");
         writeImage(diffMask, "ssdiff.mask.fits");
@@ -674 +766 @@
         kernel->xMax = SKY_SIZE;
         kernel->yMax = SKY_SIZE;
+
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernel);
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernel);
+#if 1
+//        kernel->solution1->data.F64[normIndex] += 1.0;
+        kernel->solution1->data.F64[bgIndex] = 0.0;
+#endif
+        fprintf(stderr, "Norm: %f BG: %f\n", kernel->solution1->data.F64[normIndex], kernel->solution1->data.F64[bgIndex]);
 
         pmReadout *warp = readouts->data[readouts->n - 1];