Changeset 18249

Timestamp:

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

Author:

Paul Price

Message:

Adding penalty for wideness.

Files:

• branches/pap_branch_080617/ppSub/src/ppSubArguments.c (modified) (2 diffs)
• branches/pap_branch_080617/ppSub/src/ppSubCamera.c (modified) (5 diffs)
• branches/pap_branch_080617/ppSub/src/ppSubReadout.c (modified) (5 diffs)
• branches/pap_branch_080617/ppSub/test/fake.c (modified) (4 diffs)
• trunk/ippconfig/recipes/ppSub.config (modified) (2 diffs)

branches/pap_branch_080617/ppSub/src/ppSubArguments.c

@@ -186 +186 @@
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-order", 0, "Spatial polynomial order", -1);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Kernel type (ISIS|POIS|SPAM|FRIES|GUNK|RINGS)", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-penalty", 0, "Penalty for wideness", NAN);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-widths", 0, "ISIS Gaussian FWHMs (comma-separated)", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-orders", 0, "ISIS polynomial orders (comma-separated)", NULL);
@@ -260 +261 @@
     VALUE_ARG_RECIPE_FLOAT("-rej",        "REJ",             F32);
     VALUE_ARG_RECIPE_FLOAT("-badfrac",    "BADFRAC",         F32);
+    VALUE_ARG_RECIPE_FLOAT("-penalty",    "PENALTY",         F32);
 
     valueArgRecipeStr(arguments, recipe, "-mask-bad",   "MASK.BAD",   config->arguments);

branches/pap_branch_080617/ppSub/src/ppSubCamera.c

@@ -82 +82 @@
 
     // Output image
-    pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PPSUB.OUTPUT");
+    pmFPAfile *output = pmFPAfileDefineFromFile(config, input, 1, 1, "PPSUB.OUTPUT");
     if (!output) {
         psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT"));
@@ -94 +94 @@
 
     // Output mask
-    pmFPAfile *outMask = pmFPAfileDefineSkycell(config, output->fpa, "PPSUB.OUTPUT.MASK");
+    pmFPAfile *outMask = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.OUTPUT.MASK");
     if (!outMask) {
         psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.MASK"));
@@ -107 +107 @@
     // Output weight
     if (inputWeight && refWeight) {
-        pmFPAfile *outWeight = pmFPAfileDefineSkycell(config, output->fpa, "PPSUB.OUTPUT.WEIGHT");
+        pmFPAfile *outWeight = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.OUTPUT.WEIGHT");
         if (!outWeight) {
             psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.WEIGHT"));
@@ -119 +119 @@
     }
 
+#if 0
     if (!pmFPAAddSourceFromFormat(output->fpa, "Subtraction", output->format)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA.");
         return false;
     }
+#endif
 
     pmFPAfile *sources = pmFPAfileDefineFromArgs(&status, config, "PPSUB.SOURCES", "PPSUB.SOURCES");
@@ -153 +155 @@
 
         // Internal-ish file for getting the PSF from the matched addition
-        pmFPAfile *psf = pmFPAfileDefineSkycell(config, psphot->fpa, "PSPHOT.PSF.LOAD");
+        pmFPAfile *psf = pmFPAfileDefineOutputFromFile(config, output, "PSPHOT.PSF.LOAD");
         if (!psf) {
             psError(PS_ERR_IO, false, "Failed to build FPA from PSPHOT.PSF.LOAD");

branches/pap_branch_080617/ppSub/src/ppSubReadout.c

@@ -73 +73 @@
     int ringsOrder = psMetadataLookupS32(NULL, config->arguments, "RINGS.ORDER"); // RINGS polynomial order
     int binning = psMetadataLookupS32(NULL, config->arguments, "SPAM.BINNING"); // Binning for SPAM kernel
+    float penalty = psMetadataLookupF32(NULL, config->arguments, "PENALTY"); // Penalty for wideness
     psMaskType maskBad = pmConfigMask(psMetadataLookupStr(NULL, config->arguments, "MASK.BAD"),
                                     config); // Value to mask
@@ -131 +132 @@
 
     if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, regionSize, spacing, threshold, sources,
-                            stampsName, type, size, order, widths, orders, inner, ringsOrder,
-                            binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
+                            stampsName, type, size, order, widths, orders, inner, ringsOrder, binning,
+                            penalty, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
                             maskBlank, badFrac, mode)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
@@ -184 +185 @@
 #endif
 
-    // "Subtract" the mask and weight map
+    // Subtraction is: minuend - subtrahend
+    pmReadout *minuend = inConv;
+    pmReadout *subtrahend = refConv;
+
+    if (reverse) {
+        pmReadout *temp = subtrahend;
+        subtrahend = minuend;
+        minuend = temp;
+    }
+
+#ifdef TESTING
+    {
+        psFits *fits = psFitsOpen("minuend.fits", "w");
+        psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
+        psFitsClose(fits);
+    }
+    {
+        psFits *fits = psFitsOpen("subtrahend.fits", "w");
+        psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
+            psFitsClose(fits);
+    }
+#endif
+
     outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask);
-    if (inConv->weight && refConv->weight) {
-        outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
-    }
-    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
+    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true; // It'll be there soon
 
     // Photometry is to be performed in two stages:
@@ -198 +218 @@
     pmPSF *psf = NULL;                  // PSF for photometry
     if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
-        // We use a summed image as the basis for the PSF: this will have the maximum S/N.
-        outRO->image = (psImage*)psBinaryOp(outRO->image, inConv->image, "+", refConv->image);
+        outRO->image = psImageCopy(outRO->image, minuend->image, PS_TYPE_F32);
+        if (minuend->weight) {
+            outRO->weight = psImageCopy(outRO->weight, minuend->weight, PS_TYPE_F32);
+        }
 
         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
@@ -244 +266 @@
     }
 
-    // Do the subtraction
-    {
-        // Subtraction is: minuend - subtrahend
-        pmReadout *minuend = inConv;
-        pmReadout *subtrahend = refConv;
-
-        if (reverse) {
-            pmReadout *temp = subtrahend;
-            subtrahend = minuend;
-            minuend = temp;
-        }
-
-        outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
-
-#ifdef TESTING
-        {
-            psFits *fits = psFitsOpen("minuend.fits", "w");
-            psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        {
-            psFits *fits = psFitsOpen("subtrahend.fits", "w");
-            psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
-            psFitsClose(fits);
-        }
-#endif
-
-        pmReadoutMaskApply(outRO, maskBlank);
-    }
+    // Do the actual subtraction
+    outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
+    if (inConv->weight && refConv->weight) {
+        outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
+    }
+    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
+
+    pmReadoutMaskApply(outRO, maskBlank);
 
     psFree(inConv);

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;

trunk/ippconfig/recipes/ppSub.config

@@ -18 +18 @@
 INNER           S32     5               # Inner half-size for SPAM and FRIES kernels
 RINGS.ORDER     S32     2               # Polynomial order for RINGS kernels
+PENALTY         F32     1.0             # Penalty for wideness
 
 OPTIMUM         BOOL    FALSE           # Derive optimum parameters for ISIS and GUNK kernels
@@ -29 +30 @@
 RENORM.WIDTH    S32     500             # Size of renormalisation boxes
 
-PHOTOMETRY      BOOL    TRUE            # Perform photometry?
+PHOTOMETRY      BOOL    FALSE           # Perform photometry?
 
 ### Modifications to use when stacking data