Changeset 25626 for branches/eam_branches/20090715/psModules/src/imcombine

Timestamp:

Sep 27, 2009, 11:36:17 AM (17 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/20090715/psModules/src/imcombine

Files:

• pmPSFEnvelope.c (modified) (5 diffs)
• pmStackReject.c (modified) (6 diffs)
• pmSubtractionMatch.c (modified) (1 diff)

branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c

@@ -34 +34 @@
 
 // #define TESTING                         // Enable test output
+// #define PEAK_NORM                       // Normalise peaks?
 #define PEAK_FLUX 1.0e4                 // Peak flux for each source
 #define SKY_VALUE 0.0e0                 // Sky value for fake image
@@ -122 +123 @@
             continue;
         }
+
+        if (psTraceGetLevel("psModules.imcombine") >= 1) {
+            psString string = NULL;     // String with values
+            psStringAppend(&string, "PSF %d: ", i);
+            float x = numCols / 2.0, y = numRows / 2.0; // Coordinates of interest
+            for (int j = 4; j < psf->params->n; j++) {
+                pmTrend2D *trend = psf->params->data[j]; // Trend of interest
+                double val = pmTrend2DEval(trend, x, y);
+                double err;
+                switch (trend->mode) {
+                  case PM_TREND_POLY_ORD:
+                  case PM_TREND_POLY_CHEB:
+                    err = NAN;
+                    break;
+                  case PM_TREND_MAP:
+                    err = psImageUnbinPixel(x, y, trend->map->error, trend->map->binning);
+                    break;
+                  default:
+                    psAbort("Unsupported mode: %x", trend->mode);
+                }
+                psStringAppend(&string, "%lf %lf   ", val, err);
+            }
+            psTrace("psModules.imcombine", 1, "%s\n", string);
+            psFree(string);
+        }
+
         pmResiduals *resid = psf->residuals;// PSF residuals
         psf->residuals = NULL;
@@ -144 +171 @@
             float y = source->peak->yf + yOffset->data.S32[j]; // y coordinate of source
 
-            double flux = fakeRO->image->data.F32[(int)y][(int)x];
+#ifdef PEAK_NORM
+            // Perhaps I'm being paranoid, but specify a range to check
+            int uMax = PS_MIN(x + radius, numCols - 1), uMin = PS_MAX(x - radius, 0);
+            int vMax = PS_MIN(y + radius, numRows - 1), vMin = PS_MAX(y - radius, 0);
+
+            double flux = -INFINITY;    // Peak flux
+            for (int v = vMin; v <= vMax; v++) {
+                for (int u = uMin; u <= uMax; u++) {
+                    if (fakeRO->image->data.F32[v][u] > flux) {
+                        flux = fakeRO->image->data.F32[v][u];
+                    }
+                }
+            }
             if (!isfinite(flux) || flux < 0) {
                 continue;
             }
             float norm = PEAK_FLUX / flux; // Normalisation for source
+#endif
             psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region
             psImage *subImage = psImageSubset(fakeRO->image, region); // Subimage of fake PSF
             psImage *subEnv = psImageSubset(envelope, region); // Subimage of envelope
+#ifdef PEAK_NORM
             psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32));
+#endif
             psBinaryOp(subEnv, subEnv, "MAX", subImage);
             psFree(subImage);
@@ -345 +387 @@
     psFree(try);
 
+    if (psTraceGetLevel("psModules.imcombine") >= 1) {
+        psString string = NULL;     // String with values
+        psStringAppend(&string, "Envelope PSF: ");
+        float x = numCols / 2.0, y = numRows / 2.0; // Coordinates of interest
+        for (int j = 4; j < psf->params->n; j++) {
+            pmTrend2D *trend = psf->params->data[j]; // Trend of interest
+            double val = pmTrend2DEval(trend, x, y);
+            double err;
+            switch (trend->mode) {
+              case PM_TREND_POLY_ORD:
+              case PM_TREND_POLY_CHEB:
+                err = NAN;
+                break;
+              case PM_TREND_MAP:
+                err = psImageUnbinPixel(x, y, trend->map->error, trend->map->binning);
+                break;
+              default:
+                psAbort("Unsupported mode: %x", trend->mode);
+            }
+            psStringAppend(&string, "%lf %lf   ", val, err);
+        }
+        psTrace("psModules.imcombine", 1, "%s\n", string);
+        psFree(string);
+    }
+
 #ifdef TESTING
     {
@@ -358 +425 @@
 
         pmReadout *generated = pmReadoutAlloc(NULL); // Generated image
-        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, NULL, NULL, psf, NAN, radius,
+        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, 0, NULL, NULL, psf, NAN, radius,
                                  false, true);
         {

branches/eam_branches/20090715/psModules/src/imcombine/pmStackReject.c

@@ -43 +43 @@
     if (box > 0) {
         // Convolve a subimage, then stick it in the target
-        // XXX if (threaded) {
-        // XXX     psMutexLock(source);
-        // XXX }
         psImage *mask = psImageSubset(source, psRegionSet(xMin - box, xMax + box,
                                                           yMin - box, yMax + box)); // Mask to convolve
-        // XXX if (threaded) {
-        // XXX     psMutexUnlock(source);
-        // XXX }
         psImage *convolved = psImageConvolveMask(NULL, mask, PM_STACK_MASK_BAD, PM_STACK_MASK_CONVOLVE,
                                                  -box, box, -box, box); // Convolved mask
-        // XXX if (threaded) {
-        // XXX     psMutexLock(source);
-        // XXX }
         psFree(mask);
-        // XXX if (threaded) {
-        // XXX     psMutexUnlock(source);
-        // XXX }
 
         int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_IMAGE_MASK); // Number of bytes to copy
@@ -162 +150 @@
     pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
     inRO->image = image;
-    // XXX if (threaded) {
-    // XXX     psMutexInit(image);
-    // XXX }
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = subRegions->data[i]; // Region of interest
         pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest
-        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, stride, 0, 0, 1.0, 0.0,
+        if (!pmSubtractionConvolve(NULL, convRO, NULL, inRO, NULL, stride, 0, 0, 1.0, 0.0,
                                    region, kernels, false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
@@ -211 +196 @@
         }
     }
-    // XXX if (threaded) {
-    // XXX     psMutexDestroy(image);
-    // XXX }
     psFree(inRO);
     psImage *convolved = psMemIncrRefCounter(convRO->image);
@@ -264 +246 @@
     psImage *target = psImageRecycle(convolved, numCols, numRows, PS_TYPE_IMAGE_MASK); // Grown image
     psImageInit(target, 0);
-    // XXX if (threaded) {
-    // XXX     psMutexInit(source);
-    // XXX }
     for (int i = 0; i < subRegions->n; i++) {
         psRegion *region = subRegions->data[i]; // Subtraction region
@@ -287 +266 @@
                     psArray *args = job->args; // Job arguments
                     psArrayAdd(args, 1, target);
-                    // XXX psMutexLock(source);
                     psArrayAdd(args, 1, source);
-                    // XXX psMutexUnlock(source);
                     psArrayAdd(args, 1, kernels);
                     PS_ARRAY_ADD_SCALAR(args, numCols, PS_TYPE_S32);
@@ -332 +309 @@
         }
 
-        // XXX psMutexDestroy(source);
     }
 

branches/eam_branches/20090715/psModules/src/imcombine/pmSubtractionMatch.c

@@ -628 +628 @@
     variance = NULL;
 
-    if (!pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
+    if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
+        goto MATCH_ERROR;
+    }
+    if (conv2 && !pmSubtractionBorder(conv2->image, conv2->variance, conv2->mask, size, maskBad)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
         goto MATCH_ERROR;