Index: /branches/eam_rel9_b0/psModules/src/objects/pmObjects.c
===================================================================
--- /branches/eam_rel9_b0/psModules/src/objects/pmObjects.c	(revision 5805)
+++ /branches/eam_rel9_b0/psModules/src/objects/pmObjects.c	(revision 5806)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-07 20:01:05 $
+ *  @version $Revision: 1.3.4.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-17 20:41:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -134,5 +134,6 @@
     tmp->modelPSF = NULL;
     tmp->modelFLT = NULL;
-    tmp->type = 0;
+    tmp->type = PM_SOURCE_UNKNOWN;
+    tmp->mode = PM_SOURCE_DEFAULT;
     psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
 
@@ -749,4 +750,5 @@
     psS32 numPixels = 0;
     psF32 Sum = 0.0;
+    psF32 Var = 0.0;
     psF32 X1 = 0.0;
     psF32 Y1 = 0.0;
@@ -786,4 +788,5 @@
 
             psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+            psF32 wDiff = source->weight->data.F32[row][col];
 
             // XXX EAM : check for valid S/N in pixel
@@ -793,4 +796,5 @@
             }
 
+            Var += wDiff;
             Sum += pDiff;
             X1  += xDiff * pDiff;
@@ -805,7 +809,9 @@
         }
     }
+
+    // if we have less than (1/4) of the possible pixels, force a retry
     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < 3) || (Sum <= 0)) {
-        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
+    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
+        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
         return (false);
@@ -824,5 +830,5 @@
     y = Y1/Sum;
     if ((fabs(x) > radius) || (fabs(y) > radius)) {
-        psTrace (".psModules.pmSourceMoments", 5,
+        psTrace (".psModules.pmSourceMoments", 3,
                  "large centroid swing; invalid peak %d, %d\n",
                  source->peak->x, source->peak->y);
@@ -837,4 +843,5 @@
     source->moments->Sxy = XY/Sum - x*y;
     source->moments->Sum = Sum;
+    source->moments->SN  = Sum / sqrt(Var);
     source->moments->Peak = peakPixel;
     source->moments->nPixels = numPixels;
@@ -1064,5 +1071,5 @@
 XXX: How can this function ever return FALSE?
  
-XXX EAM : add the saturated mask value to metadata 
+EAM: I moved S/N calculation to pmSourceMoments, using weight image
 *****************************************************************************/
 
@@ -1079,6 +1086,8 @@
     int Ncr      = 0;
     int Nsatstar = 0;
-    psRegion allArray = psRegionSet (0, 0, 0, 0);
-
+    // psRegion allArray = psRegionSet (0, 0, 0, 0);
+    psRegion inner;
+
+    // report stats on S/N values for star-like objects
     psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
     starsn->n = 0;
@@ -1086,8 +1095,5 @@
     // check return status value (do these exist?)
     bool status;
-    psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
-    psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
     psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
-    // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
 
     // XXX allow clump size to be scaled relative to sigmas?
@@ -1102,14 +1108,7 @@
         psF32 sigY = tmpSrc->moments->Sy;
 
-        // calculate and save signal-to-noise estimates
-        psF32 S  = tmpSrc->moments->Sum;
-        psF32 A  = 4 * M_PI * sigX * sigY;
-        psF32 B  = tmpSrc->moments->Sky;
-        psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
-        psF32 SN = (S * sqrt(GAIN) / RT);
-        tmpSrc->moments->SN = SN;
-
         // XXX EAM : can we use the value of SATURATE if mask is NULL?
-        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
+        inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
+        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
@@ -1117,5 +1116,6 @@
         bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
         if ((Nsatpix > 1) && big) {
-            tmpSrc->type = PM_SOURCE_SATSTAR;
+            tmpSrc->type = PM_SOURCE_STAR;
+            tmpSrc->mode = PM_SOURCE_SATSTAR;
             Nsatstar ++;
             continue;
@@ -1125,4 +1125,5 @@
         if (Nsatpix > 1) {
             tmpSrc->type = PM_SOURCE_SATURATED;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Nsat ++;
             continue;
@@ -1134,4 +1135,5 @@
         if ((sigX < 0.05) || (sigY < 0.05)) {
             tmpSrc->type = PM_SOURCE_DEFECT;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Ncr ++;
             continue;
@@ -1141,4 +1143,5 @@
         if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
             tmpSrc->type = PM_SOURCE_GALAXY;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Ngal ++;
             continue;
@@ -1146,12 +1149,13 @@
 
         // the rest are probable stellar objects
-        starsn->data.F32[starsn->n] = SN;
+        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
         starsn->n ++;
         Nstar ++;
 
-        // PSF star (within 1.5 sigma of clump center
+        // PSF star (within 1.5 sigma of clump center, S/N > limit)
         psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
-        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
-            tmpSrc->type = PM_SOURCE_PSFSTAR;
+        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
+            tmpSrc->type = PM_SOURCE_STAR;
+            tmpSrc->mode = PM_SOURCE_PSFSTAR;
             Npsf ++;
             continue;
@@ -1159,5 +1163,6 @@
 
         // random type of star
-        tmpSrc->type = PM_SOURCE_OTHER;
+        tmpSrc->type = PM_SOURCE_STAR;
+        tmpSrc->mode = PM_SOURCE_DEFAULT;
     }
 
@@ -1698,4 +1703,6 @@
     }
 
+    source->mode |= PM_SOURCE_FITTED;
+
     psFree(x);
     psFree(y);
@@ -1729,5 +1736,5 @@
     psBool onPic     = true;
     psBool rc        = true;
-    psF32  Ro, ymodel;
+    // psF32  Ro, ymodel;
 
     psVector *params = model->params;
@@ -1739,5 +1746,5 @@
     // XXX EAM : I need to use the sky value to constrain the weight model
     int nParams = PSF ? params->n - 4 : params->n;
-    psF32 So = params->data.F32[0];
+    // psF32 So = params->data.F32[0];
 
     // find the number of valid pixels
@@ -1776,8 +1783,8 @@
 
                 // compare observed flux to model flux to adjust weight
-                ymodel = modelFunc (NULL, model->params, coord);
+                // ymodel = modelFunc (NULL, model->params, coord);
 
                 // this test enhances the weight based on deviation from the model flux
-                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
+                // Ro = 1.0 + sqrt (fabs (y->data.F32[tmpCnt] - ymodel) / (PS_SQR(ymodel - So) + 1.0));
 
                 // psMinimizeLMChi2_EAM takes wt = 1/dY^2
@@ -1785,5 +1792,6 @@
                     yErr->data.F32[tmpCnt] = 0.0;
                 } else {
-                    yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
+                    // yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
+                    yErr->data.F32[tmpCnt] = 1.0 / source->weight->data.F32[i][j];
                 }
                 tmpCnt++;
@@ -1865,4 +1873,6 @@
         model->status = PM_MODEL_OFFIMAGE;
     }
+
+    source->mode |= PM_SOURCE_FITTED;
 
     psFree(x);