Changeset 6448 for branches/rel10_ifa/psModules/src/objects/pmObjects.c

Timestamp:

Feb 17, 2006, 7:13:42 AM (20 years ago)

Author:

magnier

Message:

bulk merge of eam_rel9_p0 onto this branch

File:

• branches/rel10_ifa/psModules/src/objects/pmObjects.c (modified) (51 diffs)

branches/rel10_ifa/psModules/src/objects/pmObjects.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-06 22:17:54 $
+ *  @version $Revision: 1.9.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 17:13:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -25 +25 @@
                     psS32 y,
                     psF32 counts,
-                    pmPeakType class)
+                    pmPeakType type)
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -32 +32 @@
     tmp->y = y;
     tmp->counts = counts;
-    tmp->class = class;
+    tmp->type = type;
 
     psTrace(__func__, 3, "---- %s() end ----\n", __func__);
@@ -77 +77 @@
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
-    psFree(tmp->modelFLT);
+    psFree(tmp->modelEXT);
+    psFree(tmp->blends);
     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
 }
@@ -85 +86 @@
 psVector containing the specified row of data from the psImage.
  
-XXX: Is there a better way to do this?
+XXX: Is there a better way to do this?  
 XXX EAM: does this really need to alloc a new vector???
 *****************************************************************************/
@@ -211 +212 @@
     // XXX EAM : i changed this to match pmModelEval above, but see
     // XXX EAM   the note below in pmSourceContour
-    psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
+    psF32 oldValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow);
     if (oldValue == level) {
         psTrace(__func__, 4, "---- %s() end ----\n", __func__);
@@ -233 +234 @@
 
     while (subCol != lastColumn) {
-        psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
+        psF32 newValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow);
         if (oldValue == level) {
             psTrace(__func__, 4, "---- %s() end ----\n", __func__);
@@ -310 +311 @@
     tmp->mask = NULL;
     tmp->moments = NULL;
+    tmp->blends = NULL;
     tmp->modelPSF = NULL;
-    tmp->modelFLT = NULL;
-    tmp->type = 0;
+    tmp->modelEXT = NULL;
+    tmp->type = PM_SOURCE_UNKNOWN;
+    tmp->mode = PM_SOURCE_DEFAULT;
     psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
 
@@ -436 +439 @@
 XXX: In the output psArray elements, should we use the image row/column offsets?
      Currently, we do not.
+XXX EAM : this function needs to return peaks in *parent* coords
  
 XXX: Merge with CVS 1.20.  This had the proper code for images with a single
 row or column.
+ 
 *****************************************************************************/
 psArray *pmFindImagePeaks(const psImage *image,
@@ -456 +461 @@
     psArray *list = NULL;
 
+    psU32 col0 = image->col0;
+    psU32 row0 = image->row0;
+
     //
     // Find peaks in row 0 only.
@@ -462 +470 @@
     tmpRow = getRowVectorFromImage((psImage *) image, row);
     psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
+    // pmFindVectorPeaks returns coords in the vector, not corrected for col0
 
     for (psU32 i = 0 ; i < row1->n ; i++ ) {
@@ -474 +483 @@
 
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
                 }
             }
@@ -484 +493 @@
                     (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
                 }
             }
@@ -493 +502 @@
                     (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
                 }
             }
@@ -532 +541 @@
                         (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                     myType = PM_PEAK_EDGE;
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
                 }
             } else if (col < (image->numCols - 1)) {
@@ -567 +576 @@
                         }
 
-                        list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
+                        list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
                     }
                 }
@@ -579 +588 @@
                         (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
                     myType = PM_PEAK_EDGE;
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
                 }
             } else {
@@ -603 +612 @@
                     (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
                 }
             }
@@ -613 +622 @@
                     (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
                 }
             }
@@ -622 +631 @@
                     (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
                 }
             }
@@ -728 +737 @@
 XXX: Sync with IfA on whether the peak x/y coords are data structure coords,
      or they use the image row/column offsets.
+XXX  EAM : peak->x,y uses parent coordinates
  
 XXX: Should we simply set pmSource->peak = peak?  If so, should we increase
@@ -776 +786 @@
         return(false);
     }
-    source->moments = pmMomentsAlloc();
+    if (source->moments == NULL) {
+        source->moments = pmMomentsAlloc();
+    }
     source->moments->Sky = (psF32) tmpF64;
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return (true);
+}
+
+// A complementary function to pmSourceLocalSky: calculate the local median variance
+bool pmSourceLocalSkyVariance(
+    pmSource *source,
+    psStatsOptions statsOptions,
+    psF32 Radius)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->weight, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_INT_POSITIVE(Radius, false);
+    PS_ASSERT_INT_NONNEGATIVE(Radius, false);
+
+    psImage *image = source->weight;
+    psImage *mask  = source->mask;
+    pmPeak *peak  = source->peak;
+    psRegion srcRegion;
+
+    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
+    srcRegion = psRegionForImage(mask, srcRegion);
+
+    psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);
+    psStats *myStats = psStatsAlloc(statsOptions);
+    myStats = psImageStats(myStats, image, mask, 0xff);
+    psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);
+
+    psF64 tmpF64;
+    p_psGetStatValue(myStats, &tmpF64);
+    psFree(myStats);
+
+    if (isnan(tmpF64)) {
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
+    }
+    if (source->moments == NULL) {
+        source->moments = pmMomentsAlloc();
+    }
+    source->moments->dSky = (psF32) tmpF64;
     psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     return (true);
@@ -830 +885 @@
     psS32 numPixels = 0;
     psF32 Sum = 0.0;
+    psF32 Var = 0.0;
     psF32 X1 = 0.0;
     psF32 Y1 = 0.0;
@@ -841 +897 @@
     psF32 xPeak = source->peak->x;
     psF32 yPeak = source->peak->y;
+    psF32 xOff = source->pixels->col0 - source->peak->x;
+    psF32 yOff = source->pixels->row0 - source->peak->y;
 
     // XXX why do I get different results for these two methods of finding Sx?
@@ -851 +909 @@
     // XXX EAM : mask == 0 is valid
 
+    psF32 **vPix = source->pixels->data.F32;
+    psF32 **vWgt = source->weight->data.F32;
+    psU8  **vMsk = source->mask->data.U8;
+
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
+            if ((source->mask != NULL) && (vMsk[row][col])) {
                 continue;
             }
 
-            psF32 xDiff = col + source->pixels->col0 - xPeak;
-            psF32 yDiff = row + source->pixels->row0 - yPeak;
+            // psF32 xDiff = col + source->pixels->col0 - xPeak;
+            // psF32 yDiff = row + source->pixels->row0 - yPeak;
+
+            psF32 xDiff = col + xOff;
+            psF32 yDiff = row + yOff;
 
             // XXX EAM : calculate xDiff, yDiff up front;
@@ -866 +931 @@
             }
 
-            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+            psF32 pDiff = vPix[row][col] - sky;
+            psF32 wDiff = vWgt[row][col];
 
             // XXX EAM : check for valid S/N in pixel
             // XXX EAM : should this limit be user-defined?
-            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
+            if (PS_SQR(pDiff) < wDiff) {
                 continue;
             }
 
+            Var += wDiff;
             Sum += pDiff;
-            X1  += xDiff * pDiff;
-            Y1  += yDiff * pDiff;
-            XY  += xDiff * yDiff * pDiff;
-
-            X2  += PS_SQR(xDiff) * pDiff;
-            Y2  += PS_SQR(yDiff) * pDiff;
-
-            peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
+
+            psF32 xWght = xDiff * pDiff;
+            psF32 yWght = yDiff * pDiff;
+
+            X1  += xWght;
+            Y1  += yWght;
+            XY  += xDiff * yWght;
+
+            X2  += xDiff * xWght;
+            Y2  += yDiff * yWght;
+
+            peakPixel = PS_MAX (vPix[row][col], peakPixel);
             numPixels++;
         }
     }
+
+    // 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);
@@ -905 +978 @@
     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);
@@ -918 +991 @@
     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;
@@ -976 +1050 @@
 
 /******************************************************************************
-pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
-sigma-x, sigma-y plane. return 0,0 clump in case of error.
+    pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 
+    sigma-x, sigma-y plane. return 0,0 clump in case of error. 
 *****************************************************************************/
 
@@ -1138 +1212 @@
 
 /******************************************************************************
-pmSourceRoughClass(source, metadata): make a guess at the source
-classification.
- 
-XXX: push the clump info into the metadata?
- 
-XXX: How can this function ever return FALSE?
- 
-XXX EAM : add the saturated mask value to metadata
+    pmSourceRoughClass(source, metadata): make a guess at the source
+    classification.
+     
+    XXX: push the clump info into the metadata?
+     
+    XXX: How can this function ever return FALSE?
+     
+    EAM: I moved S/N calculation to pmSourceMoments, using weight image
 *****************************************************************************/
 
@@ -1155 +1229 @@
 
     int Nsat     = 0;
-    int Ngal     = 0;
+    int Next     = 0;
     int Nstar    = 0;
     int Npsf     = 0;
     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;
@@ -1167 +1243 @@
     // 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?
@@ -1178 +1251 @@
         pmSource *tmpSrc = (pmSource *) sources->data[i];
 
-        tmpSrc->peak->class = 0;
+        tmpSrc->peak->type = 0;
 
         psF32 sigX = tmpSrc->moments->Sx;
         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);
+        inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2);
+        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
         // 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_SATSTAR;
+            tmpSrc->type = PM_SOURCE_STAR;
+            tmpSrc->mode = PM_SOURCE_SATSTAR;
             Nsatstar ++;
             continue;
@@ -1206 +1275 @@
         if (Nsatpix > 1) {
             tmpSrc->type = PM_SOURCE_SATURATED;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Nsat ++;
             continue;
@@ -1215 +1285 @@
         if ((sigX < 0.05) || (sigY < 0.05)) {
             tmpSrc->type = PM_SOURCE_DEFECT;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Ncr ++;
             continue;
         }
 
-        // likely unsaturated galaxy (too large to be stellar)
+        // likely unsaturated extended source (too large to be stellar)
         if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
-            tmpSrc->type = PM_SOURCE_GALAXY;
-            Ngal ++;
+            tmpSrc->type = PM_SOURCE_EXTENDED;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
+            Next ++;
             continue;
         }
 
         // 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;
@@ -1240 +1313 @@
 
         // random type of star
-        tmpSrc->type = PM_SOURCE_OTHER;
+        tmpSrc->type = PM_SOURCE_STAR;
+        tmpSrc->mode = PM_SOURCE_DEFAULT;
     }
 
@@ -1254 +1328 @@
     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
     psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:     %3d\n", Ngal);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Next:     %3d\n", Next);
     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
@@ -1264 +1338 @@
 
 /** pmSourceDefinePixels()
- *
+ * 
  * Define psImage subarrays for the source located at coordinates x,y on the
  * image set defined by readout. The pixels defined by this operation consist of
@@ -1276 +1350 @@
  * example). This function should be used to define a region of interest around a
  * source, including both source and sky pixels.
- *
+ * 
  * XXX: must code this.
- *
+ * 
  */
 bool pmSourceDefinePixels(
@@ -1294 +1368 @@
 
 /******************************************************************************
-pmSourceSetPixelsCircle(source, image, radius)
- 
-XXX: This was replaced by DefinePixels in SDRS.  Remove it.
+    pmSourceSetPixelsCircle(source, image, radius)
+     
+    XXX: This was replaced by DefinePixels in SDRS.  Remove it.
 *****************************************************************************/
 bool pmSourceSetPixelsCircle(pmSource *source,
@@ -1369 +1443 @@
 
 /******************************************************************************
-pmSourceModelGuess(source, model): This function allocates a new
-pmModel structure based on the given modelType specified in the argument list.
-The corresponding pmModelGuess function is returned, and used to
-supply the values of the params array in the pmModel structure.
- 
-XXX: Many parameters are based on the src->moments structure, which is in
-image, not subImage coords.  Therefore, the calls to the model evaluation
-functions will be in image, not subImage coords.  Remember this.
+    pmSourceModelGuess(source, model): This function allocates a new
+    pmModel structure based on the given modelType specified in the argument list.  
+    The corresponding pmModelGuess function is returned, and used to 
+    supply the values of the params array in the pmModel structure.  
+     
+    XXX: Many parameters are based on the src->moments structure, which is in
+    image, not subImage coords.  Therefore, the calls to the model evaluation
+    functions will be in image, not subImage coords.  Remember this.
 *****************************************************************************/
 pmModel *pmSourceModelGuess(pmSource *source,
@@ -1394 +1468 @@
 
 /******************************************************************************
-evalModel(source, level, row): a private function which evaluates the
-source->modelPSF function at the specified coords.  The coords are subImage, not
-image coords.
- 
-NOTE: The coords are in subImage source->pixel coords, not image coords.
- 
-XXX: reverse order of row,col args?
- 
-XXX: rename all coords in this file such that their name defines whether
-the coords is in subImage or image space.
- 
-XXX: This should probably be a public pmModules function.
- 
-XXX: Use static vectors for x.
- 
-XXX: Figure out if it's (row, col) or (col, row) for the model functions.
- 
-XXX: For a while, the first psVectorAlloc() was generating a seg fault during
-testing.  Try to reproduce that and debug.
+    evalModel(source, level, row): a private function which evaluates the
+    source->modelPSF function at the specified coords.  The coords are subImage, not
+    image coords.
+     
+    NOTE: The coords are in subImage source->pixel coords, not image coords.
+     
+    XXX: reverse order of row,col args?
+     
+    XXX: rename all coords in this file such that their name defines whether
+    the coords is in subImage or image space.
+     
+    XXX: This should probably be a public pmModules function.
+     
+    XXX: Use static vectors for x.
+     
+    XXX: Figure out if it's (row, col) or (col, row) for the model functions.
+     
+    XXX: For a while, the first psVectorAlloc() was generating a seg fault during
+    testing.  Try to reproduce that and debug.
 *****************************************************************************/
 
@@ -1441 +1515 @@
 
 /******************************************************************************
-pmSourceContour(src, img, level, mode): For an input subImage, and model, this
-routine returns a psArray of coordinates that evaluate to the specified level.
- 
-XXX: Probably should remove the "image" argument.
-XXX: What type should the output coordinate vectors consist of?  col,row?
-XXX: Why a pmArray output?
-XXX: doex x,y correspond with col,row or row/col?
-XXX: What is mode?
-XXX: The top, bottom of the contour is not correctly determined.
-XXX EAM : this function is using the model for the contour, but it should
-          be using only the image counts
+    pmSourceContour(src, img, level, mode): For an input subImage, and model, this
+    routine returns a psArray of coordinates that evaluate to the specified level.
+     
+    XXX: Probably should remove the "image" argument.
+    XXX: What type should the output coordinate vectors consist of?  col,row?
+    XXX: Why a pmArray output?
+    XXX: doex x,y correspond with col,row or row/col?
+    XXX: What is mode?
+    XXX: The top, bottom of the contour is not correctly determined.
+    XXX EAM : this function is using the model for the contour, but it should
+              be using only the image counts
 *****************************************************************************/
 psArray *pmSourceContour(pmSource *source,
@@ -1464 +1538 @@
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->modelFLT, false);
-    // XXX EAM : what is the purpose of modelPSF/modelFLT?
+    PS_ASSERT_PTR_NON_NULL(source->modelEXT, false);
+    // XXX EAM : what is the purpose of modelPSF/modelEXT?
 
     //
@@ -1558 +1632 @@
 }
 
-// XXX EAM : these are better starting values, but should be available from metadata?
-#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)
+// save a static values so they may be set externally
+static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
+static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
+
+bool pmSourceFitModelInit (float nIter, float tol)
+{
+
+    PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter;
+    PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
+    return true;
+}
+
+bool pmSourceFitModel (pmSource *source,
+                       pmModel *model,
+                       const bool PSF)
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -1579 +1655 @@
     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;
@@ -1592 +1669 @@
     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) {
+            // skip masked points
+            if (source->mask->data.U8[i][j]) {
+                continue;
+            }
+            // skip zero-weight points
+            if (source->weight->data.F32[i][j] == 0) {
+                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
+            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;
+
+    // XXX EAM : the new minimization API supplies the constraints as a struct
     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
+    psMinConstrain *constrain = psMinConstrainAlloc();
+
+    // PSF model only fits first 4 parameters, EXT model fits all
     if (PSF) {
         paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
@@ -1648 +1731 @@
         }
     }
-
-    // XXX EAM : covar must be F64?
+    constrain->paramMask = paramMask;
+
+    // Set the parameter range checks
+    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
+    modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
+
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
 
     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
 
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramMask = paramMask;
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,
-                                 x, y, yErr, modelFunc);
-    psFree(constrain);
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc);
     for (int i = 0; i < dparams->n; i++) {
         if ((paramMask != NULL) && paramMask->data.U8[i])
@@ -1673 +1756 @@
     if (paramMask != NULL) {
         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
-        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
+        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
         for (int i = 0; i < dparams->n; i++) {
             if (!paramMask->data.U8[i])
@@ -1698 +1781 @@
     }
 
+    source->mode |= PM_SOURCE_FITTED;
+
     psFree(x);
     psFree(y);
@@ -1703 +1788 @@
     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 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
-
-                // 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);
-                }
-                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
-    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-    psVector *beta_lim = NULL;
-    psVector *params_min = NULL;
-    psVector *params_max = NULL;
-
-    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
-    //           in the SDRS, I define a new psMinimization which will take these in
-    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
-    modelLimits (&beta_lim, &params_min, &params_max);
-    for (int i = 0; i < params->n; i++) {
-        covar->data.F64[0][i] = beta_lim->data.F32[i];
-        covar->data.F64[1][i] = params_min->data.F32[i];
-        covar->data.F64[2][i] = params_max->data.F32[i];
-    }
-
-    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramMask = paramMask;
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,
-                                 x, y, yErr, modelFunc);
+    psFree(constrain->paramMask);
+    psFree(constrain->paramMin);
+    psFree(constrain->paramMax);
+    psFree(constrain->paramDelta);
     psFree(constrain);
-    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]);
-    }
-
-    // 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
-    model->chisq = myMin->value;
-    model->nIter = myMin->iter;
-    model->nDOF  = y->n - nParams;
-
-    // XXX EAM 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];
-        }
-    }
-
-    // set the model success or failure status
-    if (!fitStatus) {
-        model->status = PM_MODEL_NONCONVERGE;
-    } else {
-        model->status = PM_MODEL_SUCCESS;
-    }
-
-    // XXX 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;
-    }
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(paramMask);
 
     rc = (onPic && fitStatus);