Changeset 20348

Timestamp:

Oct 23, 2008, 2:16:09 PM (18 years ago)

Author:

Paul Price

Message:

Revert to using old method of finding warp PSF: old method doesn't work properly yet.

File:

• trunk/pswarp/src/pswarpLoop.c (modified) (4 diffs)

trunk/pswarp/src/pswarpLoop.c

@@ -3 +3 @@
 
 #define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
+#define PSPHOT_FIND_PSF 0               // Use psphot's findPSF function?
+//#define TESTING 1                       // Testing output?
 
 // Lists of file rules for the detectors and skycells, and an independent list of everything else
@@ -8 +10 @@
 static char *detectorFiles[] = { "PSWARP.INPUT", "PSWARP.MASK", "PSWARP.WEIGHT", NULL };
 static char *skycellFiles[] = { "PSWARP.OUTPUT", "PSWARP.OUTPUT.MASK", "PSWARP.OUTPUT.WEIGHT", NULL };
-static char *photFiles[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT", "PSPHOT.INPUT.CMF", "PSPHOT.RESID",
+static char *photFiles[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT",
+#if PSPHOT_FIND_PSF
+                             "PSPHOT.INPUT.CMF",
+#endif
+                             "PSPHOT.RESID",
                              "PSPHOT.BACKMDL", "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
                              "PSPHOT.PSF.SAVE", "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL",
@@ -350 +356 @@
                           "Don't allow psphot to fit a fake PSF!", false);
 
-#if 0
+#if !PSPHOT_FIND_PSF
         // Get rid of the transformed sources so that we can measure the PSF
         psMetadataRemoveKey(output->analysis, "PSPHOT.SOURCES");
@@ -360 +366 @@
         pmFPAview *view = pmFPAviewAlloc(0); // View into skycell
         view->chip = view->cell = view->readout = 0;
+
+#if PSPHOT_FIND_PSF
+        // Need to adjust the weight --- the main operation in psphotReadoutFindPSF is not determining the
+        // signficance of sources, but looking at the significance of individual pixels.
+        // We can adjust the weight directly since this is a deep copy
+        pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa);
+        float vf = psMetadataLookupF32(NULL, photRO->parent->concepts, "CELL.VARFACTOR"); // Variance factor
+        psBinaryOp(photRO->weight, photRO->weight, "*", psScalarAlloc(vf, PS_TYPE_F32));
+
         if (!psphotReadoutFindPSF(config, view)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to determine PSF for warped image.");
             return false;
         }
+#else
+        if (!psphotReadout(config, view)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine PSF for warped image.");
+            return false;
+        }
+#endif
+
+#ifdef TESTING
+        {
+            #define PSF_SIZE 20         // Half-size of PSF
+            #define PSF_FLUX 10000      // Central flux for PSF
+            pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);
+            pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF");
+            psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32);
+            psImageInit(image, 0);
+            pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX);
+            pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0);
+            psFree(model);
+            psFits *fits = psFitsOpen("psf.fits", "w");
+            psFitsWriteImage(fits, NULL, image, 0, NULL);
+            psFitsClose(fits);
+            psFree(image);
+        }
+#endif
+
         psFree(view);
     }