Index: /branches/eam_rel9_b0/psModules/src/objects/pmObjects.c
===================================================================
--- /branches/eam_rel9_b0/psModules/src/objects/pmObjects.c	(revision 5828)
+++ /branches/eam_rel9_b0/psModules/src/objects/pmObjects.c	(revision 5829)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.3.4.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 20:41:29 $
+ *  @version $Revision: 1.3.4.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-23 02:23:00 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -132,4 +132,5 @@
     tmp->mask = NULL;
     tmp->moments = NULL;
+    tmp->blends = NULL;
     tmp->modelPSF = NULL;
     tmp->modelFLT = NULL;
@@ -1115,4 +1116,5 @@
         // Nsigma should be user-configured parameter
         bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
+        big = true;
         if ((Nsatpix > 1) && big) {
             tmpSrc->type = PM_SOURCE_STAR;
@@ -1571,14 +1573,7 @@
 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
-/******************************************************************************
-pmSourceFitModel(source, model): must create the appropiate arguments to the
-LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
- 
-XXX: should there be a mask value?
-XXX EAM : fit the specified model (not necessarily the one in source)
-*****************************************************************************/
-bool pmSourceFitModel_v5(pmSource *source,
-                         pmModel *model,
-                         const bool PSF)
+bool pmSourceFitModel (pmSource *source,
+                       pmModel *model,
+                       const bool PSF)
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -1589,11 +1584,12 @@
     PS_ASSERT_PTR_NON_NULL(source->mask, false);
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
+
+    // XXX EAM : is it necessary for the mask & weight to exist?  the
+    //           tests below could be conditions (!NULL)
+
     psBool fitStatus = true;
     psBool onPic     = true;
     psBool rc        = true;
 
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
     psVector *params = model->params;
     psVector *dparams = model->dparams;
@@ -1602,46 +1598,48 @@
     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
 
-    int nParams = PSF ? params->n - 4 : params->n;
-
-    // find the number of valid pixels
-    // XXX EAM : this loop and the loop below could just be one pass
-    //           using the psArrayAdd and psVectorExtend functions
-    psS32 count = 0;
+    int nParams = PSF ? 4 : params->n;
+
+    // maximum number of valid pixels
+    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAlloc(nPix);
+    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
+
+    nPix = 0;
     for (psS32 i = 0; i < source->pixels->numRows; i++) {
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                count++;
-            }
-        }
-    }
-    if (count <  nParams + 1) {
+            if (source->mask->data.U8[i][j]) {
+                continue;
+            }
+            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+            // Convert i/j to image space:
+            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+            x->data[nPix] = (psPtr *) coord;
+            y->data.F32[nPix] = source->pixels->data.F32[i][j];
+
+            // psMinimizeLMChi2 takes wt = 1/dY^2
+            if (source->weight->data.F32[i][j] == 0) {
+                continue;
+            }
+            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
+            nPix++;
+        }
+    }
+    if (nPix <  nParams + 1) {
         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
         model->status = PM_MODEL_BADARGS;
+        psFree (x);
+        psFree (y);
+        psFree (yErr);
         return(false);
     }
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAlloc(count);
-    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
-    psS32 tmpCnt = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-                // XXX: Convert i/j to image space:
-                // XXX EAM: coord order is (x,y) == (col,row)
-                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-                x->data[tmpCnt] = (psPtr *) coord;
-                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-                yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
-                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
-                //           the minimization function calculates sq()
-                tmpCnt++;
-            }
-        }
-    }
+    x->n = nPix;
+    y->n = nPix;
+    yErr->n = nPix;
 
     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
@@ -1659,161 +1657,5 @@
     }
 
-    // XXX EAM : covar must be F64?
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
-
-    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
-    for (int i = 0; i < dparams->n; i++) {
-        if ((paramMask != NULL) && paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
-    }
-
-    // save the resulting chisq, nDOF, nIter
-    model->chisq = myMin->value;
-    model->nIter = myMin->iter;
-    model->nDOF  = y->n - nParams;
-
-    // get the Gauss-Newton distance for fixed model parameters
-    if (paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
-        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
-        for (int i = 0; i < dparams->n; i++) {
-            if (!paramMask->data.U8[i])
-                continue;
-            dparams->data.F32[i] = delta->data.F64[i];
-        }
-        psFree (delta);
-    }
-
-    // set the model success or failure status
-    if (!fitStatus) {
-        model->status = PM_MODEL_NONCONVERGE;
-    } else {
-        model->status = PM_MODEL_SUCCESS;
-    }
-
-    // models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-    if (!onPic) {
-        model->status = PM_MODEL_OFFIMAGE;
-    }
-
-    source->mode |= PM_SOURCE_FITTED;
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(paramMask);
-
-    rc = (onPic && fitStatus);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
-}
-
-// XXX EAM : new version with parameter range limits and weight enhancement
-bool pmSourceFitModel (pmSource *source,
-                       pmModel *model,
-                       const bool PSF)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
-    PS_ASSERT_PTR_NON_NULL(source->weight, false);
-
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-    // psF32  Ro, ymodel;
-
-    psVector *params = model->params;
-    psVector *dparams = model->dparams;
-    psVector *paramMask = NULL;
-
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-
-    // 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];
-
-    // find the number of valid pixels
-    // XXX EAM : this loop and the loop below could just be one pass
-    //           using the psArrayAdd and psVectorExtend functions
-    psS32 count = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                count++;
-            }
-        }
-    }
-    if (count <  nParams + 1) {
-        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
-        model->status = PM_MODEL_BADARGS;
-        return(false);
-    }
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAlloc(count);
-    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
-    psS32 tmpCnt = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-                // XXX: Convert i/j to image space:
-                // XXX EAM: coord order is (x,y) == (col,row)
-                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-                x->data[tmpCnt] = (psPtr *) coord;
-                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-
-                // compare observed flux to model flux to adjust weight
-                // ymodel = modelFunc (NULL, model->params, coord);
-
-                // this test enhances the weight based on deviation from the model flux
-                // Ro = 1.0 + sqrt (fabs (y->data.F32[tmpCnt] - ymodel) / (PS_SQR(ymodel - So) + 1.0));
-
-                // psMinimizeLMChi2_EAM takes wt = 1/dY^2
-                if (source->weight->data.F32[i][j] == 0) {
-                    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];
-                }
-                tmpCnt++;
-            }
-        }
-    }
-
-    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
-                            PM_SOURCE_FIT_MODEL_TOLERANCE);
-
-    // PSF model only fits first 4 parameters, FLT model fits all
-    if (PSF) {
-        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-        for (int i = 0; i < 4; i++) {
-            paramMask->data.U8[i] = 0;
-        }
-        for (int i = 4; i < paramMask->n; i++) {
-            paramMask->data.U8[i] = 1;
-        }
-    }
-
-    // XXX EAM : I've added three types of parameter range checks
-    // XXX EAM : this requires my new psMinimization functions
+    // Set the parameter range checks
     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
     psVector *beta_lim = NULL;
@@ -1839,13 +1681,10 @@
     }
 
-    // XXX EAM: we need to do something (give an error?) if rc is false
-    // XXX EAM: psMinimizeLMChi2 does not check convergence
-
-    // XXX EAM: save the resulting chisq, nDOF, nIter
+    // save the resulting chisq, nDOF, nIter
     model->chisq = myMin->value;
     model->nIter = myMin->iter;
     model->nDOF  = y->n - nParams;
 
-    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    // get the Gauss-Newton distance for fixed model parameters
     if (paramMask != NULL) {
         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
@@ -1865,5 +1704,5 @@
     }
 
-    // XXX models can go insane: reject these
+    // models can go insane: reject these
     onPic &= (params->data.F32[2] >= source->pixels->col0);
     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);