Changeset 15814 for trunk/psModules/src/imcombine/pmPSFEnvelope.c

Timestamp:

Dec 13, 2007, 12:09:35 PM (18 years ago)

Author:

Paul Price

Message:

Checking in version that should allow ppStack to compile.

File:

• trunk/psModules/src/imcombine/pmPSFEnvelope.c (modified) (10 diffs)

trunk/psModules/src/imcombine/pmPSFEnvelope.c

@@ -72 +72 @@
     float xOrigSpacing = (float)(numCols - 2 * radius) / (float)(instances - 1); // Spacing between instances
     float yOrigSpacing = (float)(numRows - 2 * radius) / (float)(instances - 1); // Spacing between instances
-    //    int fakeSpacing = 2 * radius + 1;   // Spacing between instances (x and y) in the fake image
-    //    int fakeSize = instances * fakeSpacing; // Size of fake image
+    int fakeSpacing = 2 * radius + 1;   // Spacing between instances (x and y) in the fake image
+    int fakeSize = instances * fakeSpacing; // Size of fake image
 
     // Generate list of fake sources (instances of the PSF)
@@ -82 +82 @@
     for (int j = 0, index = 0; j < instances; j++) {
         float yOrig = j * yOrigSpacing + radius; // Source position in original image
-        //        float yFake = j * fakeSpacing + radius; // Position in fake image
-        //        int dy = yFake - yOrig;         // Difference between fake and original position
+        float yFake = j * fakeSpacing + radius; // Position in fake image
+        int dy = yFake - yOrig;         // Difference between fake and original position
 
         for (int i = 0; i < instances; i++, index++) {
             float xOrig = i * xOrigSpacing + radius; // Source position in original image
-            //            float xFake = i * fakeSpacing + radius; // Position in fake image
-            //            int dx = xFake - xOrig;     // Difference between fake and original position
+            float xFake = i * fakeSpacing + radius; // Position in fake image
+            int dx = xFake - xOrig;     // Difference between fake and original position
 
             pmSource *fake = pmSourceAlloc(); // Fake source
-            //            fake->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE);
-            fake->peak = pmPeakAlloc(xOrig, yOrig, PEAK_FLUX, PM_PEAK_LONE);
+            fake->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE);
             fake->type = PM_SOURCE_TYPE_STAR;
 
@@ -104 +103 @@
 
             fakes->data[index] = fake;
-            xOffset->data.S32[index] = 0; // dx;
-            yOffset->data.S32[index] = 0; // dy;
+            xOffset->data.S32[index] = dx;
+            yOffset->data.S32[index] = dy;
         }
     }
 
     // Generate fake images with each PSF, and take the envelope
-    psImage *envelope = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Image with envelope of PSFs
+    psImage *envelope = psImageAlloc(fakeSize, fakeSize, PS_TYPE_F32); // Image with envelope of PSFs
     psImageInit(envelope, 0.0);
     pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
     for (int i = 0; i < inputs->n; i++) {
         pmPSF *psf = inputs->data[i];   // PSF of interest
-        if (!pmReadoutFakeFromSources(fakeRO, numCols, numRows, fakes, xOffset, yOffset, psf,
-                                      NAN, radius)) {
+        pmResiduals *resid = psf->residuals;// PSF residuals
+        psf->residuals = NULL;
+        if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, xOffset, yOffset, psf,
+                                      NAN, radius, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout.");
             psFree(envelope);
@@ -122 +123 @@
             psFree(xOffset);
             psFree(fakes);
+            psf->residuals = resid;
             return NULL;
         }
+        psf->residuals = resid;
 
         // Need to renormalise sources so they all have the same peak.  You would think they do have the same
@@ -159 +162 @@
         }
 
+
+
 #ifdef TESTING
         {
@@ -183 +188 @@
 #endif
 
-    // Reset the fake source pixels to the envelope pixels
+    // Put the fake sources onto a full-size image
     pmReadout *readout = pmReadoutAlloc(NULL); // Readout to contain envelope pixels
-    readout->image = envelope;
-    readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-    psImageInit(readout->weight, WEIGHT_VAL);
+    readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    psImageInit(readout->image, 0.0);
+    for (int i = 0; i < numFakes; i++) {
+        pmSource *source = fakes->data[i]; // Fake source
+        // Position of source on fake image
+        int xFake = source->peak->x + xOffset->data.S32[i];
+        int yFake = source->peak->y + yOffset->data.S32[i];
+        psRegion region = psRegionSet(xFake - radius, xFake + radius,
+                                      yFake - radius, yFake + radius); // PSF region
+        psImage *subImage = psImageSubset(envelope, region); // Subimage of fake PSF
+
+        // Position of source on "real" image
+        int x0 = source->peak->x - radius;
+        int y0 = source->peak->y - radius;
+
+        if (!psImageOverlaySection(readout->image, subImage, x0, y0, "=")) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to overlay PSF");
+            psFree(readout);
+            psFree(xOffset);
+            psFree(yOffset);
+            psFree(fakes);
+            return NULL;
+        }
+    }
+    readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "+", psScalarAlloc(1.0, PS_TYPE_F32));
     readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
     psImageInit(readout->mask, 0);
 
+#ifdef TESTING
+    {
+        // Write out the envelope
+        psFits *fits = psFitsOpen("psf_field_full.fits", "w");
+        psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
     for (int i = 0; i < numFakes; i++) {
         pmSource *source = fakes->data[i]; // Fake source
-        float x = source->peak->xf + xOffset->data.S32[i]; // x coordinates of source
-        float y = source->peak->yf + yOffset->data.S32[i]; // y coordinates of source
+        float x = source->peak->xf;// + xOffset->data.S32[i]; // x coordinates of source
+        float y = source->peak->yf;// + yOffset->data.S32[i]; // y coordinates of source
+
+        psFree(source->pixels);
+        psFree(source->weight);
+        psFree(source->maskView);
+        psFree(source->maskObj);
+        source->pixels = NULL;
+        source->weight = NULL;
+        source->maskView = NULL;
+        source->maskObj = NULL;
 
         if (!pmSourceDefinePixels(source, readout, x, y, (float)radius)) {
@@ -205 +250 @@
         }
 
+#if 0
+        // Need to do some magic so that pmSourceMoments behaves --- the only offsets it understands are the
+        // col0,row0 in the psImage.
+        int oldCol0 = source->pixels->col0;
+        int oldRow0 = source->pixels->row0;
+        int col0 = oldCol0 - xOffset->data.S32[i];
+        int row0 = oldRow0 - yOffset->data.S32[i];
+
+        P_PSIMAGE_SET_COL0(source->pixels, col0);
+        P_PSIMAGE_SET_ROW0(source->pixels, row0);
+        P_PSIMAGE_SET_COL0(source->weight, col0);
+        P_PSIMAGE_SET_ROW0(source->weight, row0);
+        P_PSIMAGE_SET_COL0(source->maskObj, col0);
+        P_PSIMAGE_SET_ROW0(source->maskObj, row0);
+#endif
+
         if (!pmSourceMoments(source, radius)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source.");
@@ -213 +274 @@
             return NULL;
         }
-    }
+
+#if 0
+        P_PSIMAGE_SET_COL0(source->pixels, oldCol0);
+        P_PSIMAGE_SET_ROW0(source->pixels, oldRow0);
+        P_PSIMAGE_SET_COL0(source->weight, oldCol0);
+        P_PSIMAGE_SET_ROW0(source->weight, oldRow0);
+        P_PSIMAGE_SET_COL0(source->maskObj, oldCol0);
+        P_PSIMAGE_SET_ROW0(source->maskObj, oldRow0);
+#endif
+    }
+
+    psMemCheckCorruption(stderr, true);
+
 
     pmPSFOptions *options = pmPSFOptionsAlloc(); // Options for fitting a PSF
@@ -224 +297 @@
     options->radius = radius;
 
-    options->psfTrendMode = PM_TREND_POLY_ORD;
+    options->psfTrendMode = PM_TREND_MAP;
     options->psfTrendNx = 1;
     options->psfTrendNy = 1;
@@ -261 +334 @@
 
         pmReadout *generated = pmReadoutAlloc(NULL); // Generated image
-        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, xOffset, yOffset, psf, NAN, radius);
+        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, xOffset, yOffset, psf, NAN, radius,
+                                 false);
         {
             psFits *fits = psFitsOpen("psf_field_model.fits", "w");