Changeset 18249 for branches/pap_branch_080617/ppSub/test/fake.c

Timestamp:

Jun 20, 2008, 3:28:06 PM (18 years ago)

Author:

Paul Price

Message:

Adding penalty for wideness.

File:

• branches/pap_branch_080617/ppSub/test/fake.c (modified) (4 diffs)

branches/pap_branch_080617/ppSub/test/fake.c

@@ -7 +7 @@
 
 // PSF for images
-#define X_AXIS_1 5.0                    // Length of x axis (FWHM) for image 1
-#define Y_AXIS_1 3.0                    // Length of y axis (FWHM) for image 1
-#define X_AXIS_2 3.0                    // Length of x axis (FWHM) for image 1
-#define Y_AXIS_2 5.0                    // Length of y axis (FWHM) for image 1
+#define X_AXIS_1 3.0                    // Length of x axis (FWHM) for image 1
+#define Y_AXIS_1 5.0                    // Length of y axis (FWHM) for image 1
+#define X_AXIS_2 5.0                    // Length of x axis (FWHM) for image 1
+#define Y_AXIS_2 3.0                    // Length of y axis (FWHM) for image 1
 
 // Image parameters
@@ -22 +22 @@
 #define SOURCES_DENSITY 30000.0         // Source density at nominated magnitude (per deg^2)
 #define SOURCES_ZP 25.0                 // Magnitude ZP
-#define SOURCES_GAUSSIAN 0              // Gaussian sources (boolean)?
+#define SOURCES_GAUSSIAN 1              // Gaussian sources (boolean)?
 
 // Noise properties
@@ -125 +125 @@
     }
 
+#if 1
+    psImage *var1 = (psImage*)psBinaryOp(NULL, image1, "+", psScalarAlloc(PS_SQR(NOISE_1), PS_TYPE_F32));
+    psImage *var2 = (psImage*)psBinaryOp(NULL, image2, "+", psScalarAlloc(PS_SQR(NOISE_2), PS_TYPE_F32));
+#else
+    psImage *var1 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32);
+    psImage *var2 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32);
+    psImageInit(var1, PS_SQR(NOISE_1));
+    psImageInit(var2, PS_SQR(NOISE_2));
+#endif
+
     for (int y = 0; y < NUMROWS; y++) {
         for (int x = 0; x < NUMCOLS; x++) {
-            image1->data.F32[y][x] += psRandomGaussian(rng) * NOISE_1;
-            image2->data.F32[y][x] += psRandomGaussian(rng) * NOISE_2;
+            image1->data.F32[y][x] += psRandomGaussian(rng) * sqrtf(var1->data.F32[y][x]);
+            image2->data.F32[y][x] += psRandomGaussian(rng) * sqrtf(var2->data.F32[y][x]);
         }
     }
-
+#if 1
+    var1 = (psImage*)psBinaryOp(var1, image1, "+", psScalarAlloc(PS_SQR(NOISE_1), PS_TYPE_F32));
+    var2 = (psImage*)psBinaryOp(var2, image2, "+", psScalarAlloc(PS_SQR(NOISE_2), PS_TYPE_F32));
+#else
+    psImageInit(var1, PS_SQR(NOISE_1));
+    psImageInit(var2, PS_SQR(NOISE_2));
+#endif
 
     {
@@ -144 +160 @@
     }
 
-    psFree(image1);
-    psFree(image2);
-
-    psFree(rng);
-
-    psImage *noise1 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32);
-    psImage *noise2 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32);
-    psImageInit(noise1, PS_SQR(NOISE_1));
-    psImageInit(noise2, PS_SQR(NOISE_2));
     {
         psFits *fits = psFitsOpen("variance1.fits", "w");
-        psFitsWriteImage(fits, NULL, noise1, 0, NULL);
+        psFitsWriteImage(fits, NULL, var1, 0, NULL);
         psFitsClose(fits);
     }
     {
         psFits *fits = psFitsOpen("variance2.fits", "w");
-        psFitsWriteImage(fits, NULL, noise2, 0, NULL);
+        psFitsWriteImage(fits, NULL, var2, 0, NULL);
         psFitsClose(fits);
     }
-    psFree(noise1);
-    psFree(noise2);
+
+    psFree(image1);
+    psFree(image2);
+
+    psFree(var1);
+    psFree(var2);
+
+    psFree(rng);
 
     return PS_EXIT_SUCCESS;