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

Timestamp:

Oct 7, 2008, 10:12:27 AM (18 years ago)

Author:

Paul Price

Message:

Cleaning up pmSourcePSFClump, adding trace statements.

File:

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

trunk/psModules/src/objects/pmSource.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-06 13:05:13 $
+ *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-10-07 20:12:27 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -279 +279 @@
     psTrace("psModules.objects", 5, "---- begin ----\n");
 
-    # define NPIX 10
-    # define SCALE 0.1
+// XXX This should be in the recipe
+#define SCALE 0.1                       // Scale of splane
 
     psArray *peaks  = NULL;
@@ -290 +290 @@
     PS_ASSERT_PTR_NON_NULL(recipe, errorClump);
 
-    bool status = true;
-    float PSF_CLUMP_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_SN_LIM");
+    bool status = true;                 // Status of MD lookup
+    float PSF_CLUMP_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_SN_LIM");
     if (!status) {
         PSF_CLUMP_SN_LIM = 0;
@@ -298 +298 @@
     // find the sigmaX, sigmaY clump
     {
-        psStats *stats  = NULL;
-        psImage *splane = NULL;
-        int binX, binY;
-        bool status;
-
-        psF32 SX_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_SX_MAX");
-        if (!status)
+        psF32 SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
+        if (!status) {
+            psWarning("MOMENTS_SX_MAX not set in recipe");
             SX_MAX = 10.0;
-        psF32 SY_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_SY_MAX");
-        if (!status)
+        }
+        psF32 SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
+        if (!status) {
+            psWarning("MOMENTS_SY_MAX not set in recipe");
             SY_MAX = 10.0;
-        psF32 AR_MAX = psMetadataLookupF32 (&status, recipe, "MOMENTS_AR_MAX");
-        if (!status)
+        }
+        psF32 AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX");
+        if (!status) {
+            psWarning("MOMENTS_AR_MAX not set in recipe");
             AR_MAX =  3.0;
+        }
         psF32 AR_MIN = 1.0 / AR_MAX;
 
         // construct a sigma-plane image
-        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
-        psImageInit(splane, 0);  // psImageAlloc doesn't zero the data
+        int numCols = SX_MAX / SCALE, numRows = SY_MAX / 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
+        psImageInit(splane, 0);
 
         // place the sources in the sigma-plane image (ignore 0,0 values?)
-        int nValid = 0;
-        for (psS32 i = 0 ; i < sources->n ; i++)
-        {
-            pmSource *tmpSrc = (pmSource *) sources->data[i];
-            if (tmpSrc == NULL)
-                continue;
-            if (tmpSrc->moments == NULL)
-                continue;
-            if (tmpSrc->moments->SN < PSF_CLUMP_SN_LIM)
-                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;
+        int nValid = 0;                 // Number of valid sources
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *src = sources->data[i]; // Source of interest
+            if (!src || !src->moments) {
+                continue;
+            }
+
+            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->SN < PSF_CLUMP_SN_LIM) {
+                psTrace("psModules.objects", 10, "Rejecting source from clump because of low S/N (%f)\n",
+                        src->moments->SN);
+                continue;
+            }
+
+            float Mxx = src->moments->Mxx, Myy = src->moments->Myy; // Second moments
+            float ar = Mxx / Myy;       // Radius
 
             // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
-            if (fabs(tmpSrc->moments->Mxx) < 0.05)
-                continue;
-            if (fabs(tmpSrc->moments->Myy) < 0.05)
-                continue;
-            if ((tmpSrc->moments->Mxx / tmpSrc->moments->Myy) > AR_MAX)
-                continue;
-            if ((tmpSrc->moments->Mxx / tmpSrc->moments->Myy) < AR_MIN)
-                continue;
-            if (tmpSrc->mode & PM_SOURCE_MODE_BLEND)
-                continue;
+            if (fabs(Mxx) < 0.05 || fabs(Myy < 0.05)) {
+                psTrace("psModules.objects", 10,
+                        "Rejecting source from clump because of low moments (%f,%f)\n",
+                        Mxx, Myy);
+                continue;
+            }
+            if (Mxx > SX_MAX || Myy > SY_MAX) {
+                psTrace("psModules.objects", 10,
+                        "Rejecting source from clump because of high moments (%f,%f)\n",
+                        Mxx, Myy);
+                continue;
+            }
+            if (ar > AR_MAX || ar < AR_MIN) {
+                psTrace("psModules.objects", 10, "Rejecting source from clump because of Ar (%f)\n", ar);
+                continue;
+            }
 
             // for the moment, force splane dimensions to be 10x10 image pix
-            binX = tmpSrc->moments->Mxx/SCALE;
-            if (binX < 0)
-                continue;
-            if (binX >= splane->numCols)
-                continue;
-
-            binY = tmpSrc->moments->Myy/SCALE;
-            if (binY < 0)
-                continue;
-            if (binY >= splane->numRows)
-                continue;
+            int binX = Mxx / SCALE, binY = Myy /  SCALE; // Position on splane
+            psAssert(binX >= 0 && binX < numCols && binY >= 0 && binY < numRows, "We checked it already");
 
             splane->data.F32[binY][binX] += 1.0;
-            nValid ++;
+            nValid++;
         }
 
         // find the peak in this image
-        stats = psStatsAlloc (PS_STAT_MAX);
+        psStats *stats = psStatsAlloc (PS_STAT_MAX);
         if (!psImageStats (stats, splane, NULL, 0)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
@@ -433 +442 @@
                 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 (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)
@@ -521 +530 @@
         pmSource *source = (pmSource *) sources->data[i];
 
-        if (source->peak->x < region->x0) continue;
-        if (source->peak->x > region->x1) continue;
-        if (source->peak->y < region->y0) continue;
-        if (source->peak->y > region->y1) continue;
+        if (source->peak->x < region->x0) continue;
+        if (source->peak->x > region->x1) continue;
+        if (source->peak->y < region->y0) continue;
+        if (source->peak->y > region->y1) continue;
 
         source->peak->type = 0;
@@ -721 +730 @@
                 vMsk++;
             }
-            if (isnan(*vPix)) continue;
+            if (isnan(*vPix)) continue;
 
             psF32 xDiff = col + xOff;
@@ -893 +902 @@
         psF32 **target = source->pixels->data.F32;
         if (mode & PM_MODEL_OP_NOISE) {
-            // XXX need to scale by the gain...
+            // XXX need to scale by the gain...
             target = source->weight->data.F32;
         }
@@ -967 +976 @@
         return model;
 
-        // the 'best' extended model is saved in source->modelEXT (may be a pointer to one of
-        // the elements of source->modelFits)
+        // the 'best' extended model is saved in source->modelEXT (may be a pointer to one of
+        // the elements of source->modelFits)
       case PM_SOURCE_TYPE_EXTENDED:
-        model = source->modelEXT;
+        model = source->modelEXT;
         if (!model && source->modelPSF) {
-            // XXX raise an error or warning here?
+            // XXX raise an error or warning here?
             if (isPSF) {
                 *isPSF = true;