Changeset 25754 for trunk/psModules/src/objects/pmSource.c

Timestamp:

Oct 2, 2009, 3:11:32 PM (17 years ago)

Author:

eugene

Message:

check in changes from genes development branch : extensive changes to moments calculation, psf model generation, aperture residuals

File:

• trunk/psModules/src/objects/pmSource.c (modified) (10 diffs)

trunk/psModules/src/objects/pmSource.c

@@ -282 +282 @@
 
     psArray *peaks  = NULL;
-    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0};
-    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0};
+    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0, 0, 0.0};
+    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0, 0, 0.0};
     pmPSFClump psfClump;
 
@@ -319 +319 @@
 
         // construct a sigma-plane image
-        int numCols = SX_MAX / PSF_CLUMP_GRID_SCALE, numRows = SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
+        int numCols = 1 + SX_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
+        int numRows = 1 + SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
         psTrace("psModules.objects", 10, "sigma-plane dimensions: %dx%d\n", numCols, numRows);
         psImage *splane = psImageAlloc(numCols, numRows, PS_TYPE_F32); // sigma-plane image
@@ -332 +333 @@
             }
 
-            int x = src->peak->x, y = src->peak->y; // Coordinates of peak
-            if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
-                continue;
-            }
+            if (region) {
+                int x = src->peak->x, y = src->peak->y; // Coordinates of peak
+                if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
+                    continue;
+                }
+            }
 
             if (src->mode & PM_SOURCE_MODE_BLEND) {
                 continue;
             }
+
+            if (!src->moments->nPixels) continue;
 
             if (src->moments->SN < PSF_SN_LIM) {
@@ -384 +389 @@
 
         // find the peak in this image
-        psStats *stats = psStatsAlloc (PS_STAT_MAX);
+        psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_SAMPLE_STDEV);
         if (!psImageStats (stats, splane, NULL, 0)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
@@ -393 +398 @@
         peaks = pmPeaksInImage (splane, stats[0].max / 2);
         psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2);
+
+        psfClump.nSigma = stats->sampleStdev;
 
         const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
@@ -430 +437 @@
         psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value);
 
+        // XXX store the mean sigma?
+        float meanSigma = psfClump.nSigma;
+        psfClump.nStars = clump->value;
+        psfClump.nSigma = clump->value / meanSigma;
+
         // define section window for clump
         minSx = clump->x * PSF_CLUMP_GRID_SCALE - 2.0*PSF_CLUMP_GRID_SCALE;
@@ -452 +464 @@
                 continue;
 
-            if (tmpSrc->peak->x < region->x0) continue;
-            if (tmpSrc->peak->x > region->x1) continue;
-            if (tmpSrc->peak->y < region->y0) continue;
-            if (tmpSrc->peak->y > region->y1) continue;
+            if (region) {
+                if (tmpSrc->peak->x < region->x0) continue;
+                if (tmpSrc->peak->x > region->x1) continue;
+                if (tmpSrc->peak->y < region->y0) continue;
+                if (tmpSrc->peak->y > region->y1) continue;
+            }
 
             if (tmpSrc->moments->Mxx < minSx)
@@ -919 +933 @@
     if (model == NULL) return false;  // model must be defined
 
+    bool addNoise = mode & PM_MODEL_OP_NOISE;
+
     if (source->modelFlux) {
         // add in the pixels from the modelFlux image
@@ -940 +956 @@
 
         psF32 **target = source->pixels->data.F32;
-        if (mode & PM_MODEL_OP_NOISE) {
-            // XXX need to scale by the gain...
+        if (addNoise) {
+            // when adding noise, we assume the shape and Io have been modified
             target = source->variance->data.F32;
         }
 
-        // XXX need to respect the source and model masks?
         for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
             int oy = iy + dY;
@@ -959 +974 @@
             }
         }
+        if (!addNoise) {
+            if (add) {
+                source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+            } else {
+                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+            }
+        }
         return true;
     }
 
     psImage *target = source->pixels;
-    if (mode & PM_MODEL_OP_NOISE) {
+    if (addNoise) {
         target = source->variance;
     }
 
-    if (add) {
-        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
-    } else {
-        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+    if (!addNoise) {
+        if (add) {
+            status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+        } else {
+            status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+        }
     }