Changeset 4954

Timestamp:

Sep 6, 2005, 5:33:01 PM (21 years ago)

Author:

eugene

Message:

psphot reorganization: merged in pmObjects code from psModule, modified to work with psLib 0.7.0

Location:

trunk/psphot

Files:

• Makefile (modified) (3 diffs)
• doc/notes.txt (modified) (1 diff)
• doc/psphot.txt (modified) (1 diff)
• src/fitsource.c (deleted)
• src/fs_args.c (deleted)
• src/models (added)
• src/models/pmModel_GAUSS.c (added)
• src/models/pmModel_PGAUSS.c (added)
• src/models/pmModel_QGAUSS.c (added)
• src/models/pmModel_RGAUSS.c (added)
• src/models/pmModel_SGAUSS.c (added)
• src/models/pmModel_TGAUSS.c (added)
• src/models/pmModel_WAUSS.c (added)
• src/models/pmModel_ZGAUSS.c (added)
• src/onesource.c (deleted)
• src/pmModelInit.c (added)
• src/pmObjects_EAM.c (modified) (51 diffs)
• src/pmObjects_EAM.h (modified) (3 diffs)
• src/pmPeaksSigmaLimit.c (modified) (1 diff)
• src/pmSourceUtils.c (modified) (1 diff)
• src/psEllipse.c (added)
• src/psEllipse.h (added)
• src/psImageData.c (modified) (1 diff)
• src/psLibUtils.c (modified) (15 diffs)
• src/psLine.c (modified) (1 diff)
• src/psMinimize.c (modified) (1 diff)
• src/psModulesUtils.c (modified) (4 diffs)
• src/psPolynomials.c (modified) (16 diffs)
• src/psphot.h (modified) (5 diffs)
• src/psphotApplyPSF.c (modified) (1 diff)
• src/psphotArguments.c (modified) (1 diff)
• src/psphotFitGalaxies.c (modified) (2 diffs)
• src/psphotOutput.c (modified) (2 diffs)
• src/psphotSetup.c (modified) (1 diff)
• src/psphotSourceStats.c (modified) (1 diff)
• src/pspsf.c (modified) (5 diffs)

trunk/psphot/Makefile

@@ -14 +14 @@
 
 INCS    =       $(IPSLIB)
-LIBS    =       -lpsmodule $(LPSLIB)
-CFLAGS  =       $(INCS) -std=c99
+# LIBS  =       -lpsmodule $(LPSLIB)
+LIBS    =       $(LPSLIB)
+CFLAGS  =       $(INCS) -std=c99 -Wall -Werror
 LFLAGS  =       $(LIBS) 
 
@@ -30 +31 @@
 $(SRC)/psphotMarkPSF.$(ARCH).o     \
 $(SRC)/psphotSubtractPSF.$(ARCH).o \
+$(SRC)/psphotSortBySN.$(ARCH).o    \
 $(SRC)/pspsf.$(ARCH).o             \
 $(SRC)/psLibUtils.$(ARCH).o        \
@@ -39 +41 @@
 $(SRC)/psMinimize.$(ARCH).o        \
 $(SRC)/psImageData.$(ARCH).o       \
-$(SRC)/psphotSortBySN.$(ARCH).o
+$(SRC)/psEllipse.$(ARCH).o         \
+$(SRC)/pmModelInit.$(ARCH).o       \
+$(SRC)/pmObjects_EAM.$(ARCH).o
 
 FITSOURCE = \

trunk/psphot/doc/notes.txt

@@ -1 +1 @@
 
 Notes on psphot
+
+2005.09.06
+
+  I have built a working version of PSPhot using my own copy of the
+  pmObjects.[ch] files.  I need to minimize the difference between the
+  v0.7.0, v0.5.0, and my own version of pmObjects.c
+
+  changes v0.5.0 / v0.7.0 / EAM
+  - changed all ps* to pm*
+  - fixed naming for local static functions (eg modelFree not p_psModelFree)
+  - dropped all hard-coded model names (vs v0.7.0)
+  - using current psMemSetDellocator name
+  - added noise image to pmSource
+  - converted old asserts to new asserts (PS_.._CHECK..)
+  - add return from error in pmFindImagePeaks
+  - allow isItInThisRegion to have NULL region (true)
+  - fixed subImages to correspond to region definition (end is
+    exclusive) in LocalSky and SetPixelCircle
+  - dropped old, redundant code in pmSourceLocalSky
+  - fixed pmSourceMoments use of mask (0 is valid, not 1)
+  - added pmSourceMoments validity tests (numPixels, Sum, centroid shift)
+  - added correction for Sxy (XY/Sum - x*y)
+  - added pmSourcePSFClump function (and pmPSFClump structure)
+  - allow missing sources in pmSourcePSFClump (not all peaks yield
+    sources)
+  - moved pmSourceRoughClass metadata lookups out of inner loop
+  - added saturated pixel test using pmImageCountPixelMask
+  - various redefinitions of initial classes
+  - dropped old, redundant code in pmSourceSetPixelsCircle
+  - mask image set to type psU8 (was psF32)
+  - abstract object model functions for pmSourceModelGuess
+  - converted evalModel to a public function psModelEval, modified params
+  - pmSourceFitModel operates on a specified model
+  - pmSourceFitModel had option to apply/ignore PSF information
+  - pmSourceFitModel calculates yErr
+  - pmSourceFitModel calculates the covariance matrix
+  - added paramMask to pmSourceFitModel 
+  - pmSourceFitModel_EAM uses covar to carry in beta and param
+    limits...
+  - pmSourceFitModel requires solution to stay within image
+  - pmSourceFitModel calculates and saves chisq, nIter, nDOF
+  - pmSourceFitModel uses GaussNewtonDelta for frozen params
+  - added mask to pmSourceAdd, pmSourceSub
+  - model abstractions in pmSourceAdd, pmSourceSub
+  - fixed 'center' option in pmSourceAdd/Sub
+
+  - pmSourceFitModel calculates sqrt(var), psMinimize calculates sq(sqrt(var))
+  - pmSourceContour is using the model flux, not the image flux
+  - does psImageAlloc zero the image?  if so, drop init routine
+    from pmSourcePSFClump
+  - pmSourceRoughClass uses RDNOISE, GAIN to calculate SN: should use
+    the source->noise image? part of pmSourceMoments?
+  - pmSourceMoments_EAM uses macro for checkRegion
+  - pmObjects_EAM.c uses some hardcoded mask values (LocalSky)
+  - why does LocalSky_EAM not need subImageMask->col0 = subImage->col0??
+  - psGetRowVectorFromImage should be in psLib (in psLibUtils)
+  - pmSourceSetPixelsCircle and pmSourceLocalSky probably could use
+    the psImageMaskRegion, etc functions.
+  - pmSourceSetPixelsCircle uses a hard-coded mask value
+  - 
 
 2005.09.05

trunk/psphot/doc/psphot.txt

@@ -1 +1 @@
 
-psphot is a stand-alone program which implements the Pan-STARRS IPP
-object detection and analysis functions for using with individual
-images.  
+Defined APIs:
 
+pmMoments           *pmMomentsAlloc();
+pmModel             *pmModelAlloc();
+pmSource            *pmSourceAlloc();
+psVector            *pmFindVectorPeaks()
+psArray             *pmFindImagePeaks()
+psList              *pmCullPeaks()
+pmSource            *pmSourceLocalSky()
+bool                 pmSourceMoments()
+pmPSFClump           pmSourcePSFClump()
+bool                 pmSourceRoughClass()
+bool                 pmSourceSetPixelsCircle()
+pmModel             *pmSourceModelGuess()
+psArray             *pmSourceContour()
+bool                 pmSourceFitModel()
+bool                 pmSourceAddModel()
+bool                 pmSourceSubModel()
+int                  pmModelParameterCount ();
+char                *pmModelGetType ();
+pmModelType          pmModelSetType ();
+pmModelFunc          pmModelFunc_GetFunction ();
+pmModelFlux          pmModelFlux_GetFunction ();
+pmModelRadius        pmModelRadius_GetFunction ();
+pmModelLimits        pmModelLimits_GetFunction ();
+pmModelGuessFunc     pmModelGuessFunc_GetFunction ();
+pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction ();
+pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction ();
+
+model function abstractions:
+  typedef psMinimizeLMChi2Func pmModelFunc;
+  typedef psF64 (*pmModelFlux)(const psVector *params);
+  typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
+  typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
+  typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
+  typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
+  typedef bool (*pmModelFitStatusFunc)(pmModel *model);

trunk/psphot/src/pmObjects_EAM.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-06 08:05:08 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-07 03:33:01 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -66 +66 @@
 pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
 and initialize the type member.  Initialize the params to 0.0.
-XXX EAM: changing params and dparams to psVector
-XXX EAM: simplifying code with psModelParameterCount
+XXX EAM: simplifying code with pmModelParameterCount
 *****************************************************************************/
 pmModel *pmModelAlloc(pmModelType type)
@@ -76 +75 @@
     tmp->chisq = 0.0;
     tmp->nIter = 0;
-    psS32 Nparams = psModelParameterCount (type);
+    psS32 Nparams = pmModelParameterCount (type);
     if (Nparams == 0) {
-        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
         return(NULL);
     }
@@ -167 +166 @@
     //
     if ((vector->data.F32[0] > vector->data.F32[1]) &&
-            (vector->data.F32[0] > threshold)) {
+        (vector->data.F32[0] > threshold)) {
         count++;
     }
@@ -176 +175 @@
     for (psU32 i = 1; i < n-1 ; i++) {
         if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
-                (vector->data.F32[i] > vector->data.F32[i+1]) &&
-                (vector->data.F32[i] > threshold)) {
+            (vector->data.F32[i] > vector->data.F32[i+1]) &&
+            (vector->data.F32[i] > threshold)) {
             count++;
         }
@@ -186 +185 @@
     //
     if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
-            (vector->data.F32[n-1] > threshold)) {
+        (vector->data.F32[n-1] > threshold)) {
         count++;
     }
@@ -201 +200 @@
     //
     if ((vector->data.F32[0] > vector->data.F32[1]) &&
-            (vector->data.F32[0] > threshold)) {
+        (vector->data.F32[0] > threshold)) {
         tmpVector->data.U32[count++] = 0;
     }
@@ -210 +209 @@
     for (psU32 i = 1; i < (n-1) ; i++) {
         if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
-                (vector->data.F32[i] > vector->data.F32[i+1]) &&
-                (vector->data.F32[i] > threshold)) {
+            (vector->data.F32[i] > vector->data.F32[i+1]) &&
+            (vector->data.F32[i] > threshold)) {
             tmpVector->data.U32[count++] = i;
         }
@@ -220 +219 @@
     //
     if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
-            (vector->data.F32[n-1] > threshold)) {
+        (vector->data.F32[n-1] > threshold)) {
         tmpVector->data.U32[count++] = n-1;
     }
@@ -319 +318 @@
         if (col == 0) {
             if ( (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
+                 (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
+                 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
 
                 if (image->data.F32[row][col] > threshold) {
@@ -328 +327 @@
         } else if (col < (image->numCols - 1)) {
             if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
+                 (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
+                 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
+                 (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
+                 (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);
@@ -339 +338 @@
         } else if (col == (image->numCols - 1)) {
             if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
-                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
+                 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                 (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);
@@ -375 +374 @@
                 // If col==0, then we can not read col-1 pixels
                 if ((image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
+                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
+                    (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);
@@ -385 +384 @@
                 // This is an interior pixel
                 if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                        (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
+                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                     if (image->data.F32[row][col] > threshold) {
                         if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row-1][col]) &&
-                                (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
-                                (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
-                                (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
+                            (image->data.F32[row][col] > image->data.F32[row-1][col]) &&
+                            (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
+                            (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
+                            (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
+                            (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
+                            (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                            (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
                             myType = PM_PEAK_LONE;
                         }
 
                         if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row-1][col]) ||
-                                (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||
-                                (image->data.F32[row][col] == image->data.F32[row+1][col]) ||
-                                (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {
+                            (image->data.F32[row][col] == image->data.F32[row-1][col]) ||
+                            (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||
+                            (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
+                            (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
+                            (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||
+                            (image->data.F32[row][col] == image->data.F32[row+1][col]) ||
+                            (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {
                             myType = PM_PEAK_FLAT;
                         }
@@ -421 +420 @@
                 // If col==numCols - 1, then we can not read col+1 pixels
                 if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                        (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
-                        (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
+                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
+                    (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
+                    (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);
@@ -448 +447 @@
         if (col == 0) {
             if ( (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
+                 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
+                 (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);
@@ -456 +455 @@
         } else if (col < (image->numCols - 1)) {
             if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col-1]) &&
-                    (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
+                 (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
+                 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
+                 (image->data.F32[row][col] >  image->data.F32[row][col-1]) &&
+                 (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);
@@ -467 +466 @@
         } else if (col == (image->numCols - 1)) {
             if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
-                    (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
+                 (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
+                 (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);
@@ -495 +494 @@
         (y >= valid->y0) &&
         (y <= valid->y1)) {
-      return(true);
+        return(true);
     }
 
@@ -511 +510 @@
  
 XXX: changed API to create a NEW output psArray (should change name as well)
+ 
+XXX: Do we free the psList elements of those culled peaks?
+
+XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?
 *****************************************************************************/
 psList *pmCullPeaks(psList *peaks,
@@ -517 +520 @@
 {
     PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-    //    PS_ASSERT_PTR_NON_NULL(valid, NULL);
 
     psListElem *tmpListElem = (psListElem *) peaks->head;
@@ -526 +528 @@
         pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;
         if ((tmpPeak->counts > maxValue) ||
-                ((valid != NULL) &&
-                 (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) {
+            ((valid != NULL) &&
+             (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) {
             psListRemoveData(peaks, (psPtr) tmpPeak);
         }
@@ -625 +627 @@
     // XXX EAM : I added this code to stay on the image. So did George
     // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
+    // XXX EAM : dropped the Annulus width (not needed)
+    // XXX EAM : dropped off-boundary tests (not needed)
     psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - outerRadiusS32);
     psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + outerRadiusS32 + 1);
     psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - outerRadiusS32);
     psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + outerRadiusS32 + 1);
-    // AnulusWidth == number of pixels width in the annulus.  We add one since
-    // the pixels at the inner AND outher radius are included.
-    // XXX EAM : not used : psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32);
-    // Example: assume an outer/inner radius of 20/10.  Then the subimage
-    // should have width/length of 40.  An 18-by-18 interior region will
-    // be masked.
-    //    printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32);
-    //    printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32);
-    //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
-
-    // XXX EAM : these tests should not be needed: we can never hit this error because of above
-    # if (1)
-
-        if (SubImageStartRow < 0) {
-            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                    SubImageStartRow);
-            return(NULL);
-        }
-    if (SubImageEndRow > image->numRows) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
-                SubImageEndRow);
-        return(NULL);
-    }
-    if (SubImageStartCol < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
-                SubImageStartCol);
-        return(NULL);
-    }
-    if (SubImageEndCol > image->numCols) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
-                SubImageEndCol);
-        return(NULL);
-    }
-    # endif
-
     //
     // Grab a subimage of the original image of size (2 * outerRadius).
     //
-    // XXX: Must fix for new psImageSubset
-    //    psImage *subImage = psImageSubset((psImage *) image,
-    //                                      SubImageStartCol,
-    //                                      SubImageStartRow,
-    //                                      SubImageEndCol,
-    //                                      SubImageEndRow);
-    //    printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows);
-    psRegion tmpRegion = psRegionSet(SubImageStartCol,
-                                     SubImageEndCol,
-                                     SubImageStartRow,
-                                     SubImageEndRow);
-    psImage *subImage = psImageSubset((psImage *) image, tmpRegion);
-    // XXX EAM : can merge these: psImageSubset (image, psRegionSet ());
+    // XXX EAM : merged psImageSubset & psRegionSet
+    // XXX EAM : cast for image is needed because of const (above)
+    psImage *subImage = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
+                                                         SubImageEndCol,
+                                                         SubImageStartRow,
+                                                         SubImageEndRow));
 
     psImage *subImageMask = psImageAlloc(subImage->numCols,
@@ -685 +647 @@
                                          PS_TYPE_U8);
     // XXX EAM : for consistency, mask needs col0,row0 set to match image
-    subImageMask->col0 = subImage->col0;
-    subImageMask->row0 = subImage->row0;
-
-    //
-    // Loop through the subimage mask, initialize mask to 1 (invalid pixel)
+    // XXX EAM : CONFLICT between psLib and code requirement. FIX PSLIB!!!
+    // subImageMask->col0 = subImage->col0;
+    // subImageMask->row0 = subImage->row0;
+    psAbort ("pmObjects", "must fix this error before psphot will work!!!\n");
+
+    //
+    // Loop through the subimage mask, initialize mask to 0x01 (invalid pixel)
     // XXX EAM : use PSPHOT_MASK_INVALID?
     //
     for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
         for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
-            subImageMask->data.U8[row][col] = 1;
+            subImageMask->data.U8[row][col] = 0x01;
         }
     }
@@ -702 +666 @@
     // XXX this uses a static mask value of 0
     // XXX EAM : this was wrong: it masked the wrong pixels at the edge of the image
+    // XXX EAM : this masks the pixels in the center region
+    // XXX EAM : should we be using psImageMaskRegion ???
     psS32 StartRow  = PS_MAX (0, SubImageCenterRow - subImageMask->row0 - innerRadiusS32);
     psS32 EndRow    = PS_MIN (subImageMask->numRows, SubImageCenterRow - subImageMask->row0 + innerRadiusS32 + 1);
@@ -712 +678 @@
     }
 
-
+    // XXX EAM : make this trace information??
     //    for (psS32 row = 0 ; row < subImage->numRows; row++) {
     //        for (psS32 col = 0 ; col < subImage->numCols; col++) {
@@ -800 +766 @@
 XXX: mask values?
 *****************************************************************************/
-bool pmSourceMoments(psSource *source,
+bool pmSourceMoments(pmSource *source,
                      psF32 radius)
 {
@@ -849 +815 @@
             psS32 imgColCoord = col + source->pixels->col0;
             psS32 imgRowCoord = row + source->pixels->row0;
-            if (CheckRadius(source->peak,
+            if (checkRadius(source->peak,
                             radius,
                             imgColCoord,
@@ -874 +840 @@
     // 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");
-      return (false);
+        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
+        return (false);
     }
 
@@ -890 +856 @@
     y = Y1/Sum;
     if ((fabs(x) > radius) || (fabs(y) > radius)) {
-      psTrace (".psModules.pmSourceMoments", 5, 
-               "large centroid swing; invalid peak %d, %d\n", 
-               source->peak->x, source->peak->y);
-      return (false);
+        psTrace (".psModules.pmSourceMoments", 5, 
+                 "large centroid swing; invalid peak %d, %d\n", 
+                 source->peak->x, source->peak->y);
+        return (false);
     }
 
@@ -914 +880 @@
 
     return(true);
-
-    // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it
-    # if (0)
-        //
-        // second loop: get the difference sums
-        //
-        X2 = Y2 = 0;
-    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] != 0)) {
-                psS32 imgColCoord = col + source->pixels->col0;
-                psS32 imgRowCoord = row + source->pixels->row0;
-                if (checkRadius(source->peak,
-                                radius,
-                                imgColCoord,
-                                imgRowCoord)) {
-                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
-                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
-                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
-
-                    Sum+= pDiff;
-                    X2+= PS_SQR(xDiff - x) * pDiff;
-                    Y2+= PS_SQR(yDiff - y) * pDiff;
-                }
-            }
-        }
-    }
-
-    //
-    // second moment X = sqrt (X2/Sum)
-    //
-    source->moments->Sx = (X2/Sum);
-    source->moments->Sy = (Y2/Sum);
-    return(true);
-    # endif
 }
 
@@ -990 +921 @@
 {
 
-    # define NPIX 10
-    # define SCALE 0.1
+# define NPIX 10
+# define SCALE 0.1
 
     psArray *peaks  = NULL;
@@ -997 +928 @@
     pmPSFClump psfClump = emptyClump;
 
-    //    PS_PTR_CHECK_NULL(sources, emptyClump);
-    //    PS_PTR_CHECK_NULL(metadata, emptyClump);
+    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
+    PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
 
     // find the sigmaX, sigmaY clump
@@ -1010 +941 @@
         splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
         for (int i = 0; i < splane->numRows; i++) {
-          memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
+            memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
         }
 
@@ -1139 +1070 @@
 
 XXX: How can this function ever return FALSE?
+
+XXX EAM : add the saturated mask value to metadata 
 *****************************************************************************/
 
@@ -1178 +1111 @@
         psF32 A  = 4 * M_PI * sigX * sigY;
         psF32 B  = tmpSrc->moments->Sky;
-        psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN)));
-        psF32 SN = (S * PS_SQRT_F32(GAIN) / RT);
+        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 SATIRATE if mask is NULL?
-        int Nsatpix = pmCountSatPixels (tmpSrc->mask);
+        // XXX EAM : can we use the value of SATURATE if mask is NULL?
+# define MASK_SATURATED 0x02
+        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
@@ -1253 +1187 @@
 }
 
-int pmCountSatPixels (psImage *image) 
-{
-    int Nsatpix = 0;
-  
-    for (int i = 0; i < image->numRows; i++) {
-        for (int j = 0; j < image->numCols; j++) {
-            if (image->data.U8[i][j] & 0x02) {
-                Nsatpix ++;
-            }
-        }
-    }
-    return (Nsatpix);
-}
-
 /******************************************************************************
 pmSourceSetPixelsCircle(source, image, radius)
@@ -1286 +1206 @@
     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->peak, false);
     PS_FLOAT_COMPARE(0.0, radius, false);
 
@@ -1305 +1225 @@
     psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
 
-    // XXX EAM : this should not be needed: we can never hit this error
-    # if (1)
-
-        if (SubImageStartRow < 0) {
-            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                    SubImageStartRow);
-            return(false);
-        }
-    if (SubImageEndRow > image->numRows) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
-                SubImageEndRow);
-        return(false);
-    }
-    if (SubImageStartCol < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
-                SubImageStartCol);
-        return(false);
-    }
-    if (SubImageEndCol > image->numCols) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
-                SubImageEndCol);
-        return(false);
-    }
-    # endif
-
     // XXX: Must recycle image.
     // XXX EAM: this message reflects a programming error we know about.
@@ -1338 +1233 @@
         psFree(source->pixels);
     }
-    // XXX: Must fix this.  psImageSubset() has different parameters in latest CVS.
-    //    source->pixels = psImageSubset((psImage *) image,
-    //                                   SubImageStartCol,
-    //                                   SubImageStartRow,
-    //                                   SubImageEndCol,
-    //                                   SubImageEndRow);
     source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
                                                                   SubImageStartRow,
@@ -1390 +1279 @@
 functions will be in image, not subImage coords.  Remember this.
 *****************************************************************************/
-psModel *pmSourceModelGuess(psSource *source,
-                            psModelType modelType)
-{
-    PS_PTR_CHECK_NULL(source, false);
-    PS_PTR_CHECK_NULL(source->moments, false);
-    PS_PTR_CHECK_NULL(source->peak, false);
-
-    psModel *model = psModelAlloc(modelType);
-
-    psModelGuessFunc modelGuessFunc = psModelGuessFunc_GetFunction (modelType);
+pmModel *pmSourceModelGuess(pmSource *source,
+                            pmModelType modelType)
+{
+    PS_ASSERT_PTR_NON_NULL(source->moments, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+
+    pmModel *model = pmModelAlloc(modelType);
+
+    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction (modelType);
     modelGuessFunc (model, source);
     return(model);
@@ -1425 +1313 @@
 testing.  Try to reproduce that and debug.
 *****************************************************************************/
-static psF32 evalModel(pmSource *src,
-                       psU32 row,
-                       psU32 col)
-{
-    PS_ASSERT_PTR_NON_NULL(src, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false);
-
-XXX EAM : I've made this a public function
-XXX EAM : this now uses a psModel as the input
-XXX EAM : it was using src->type to find the model, not model->type
-*****************************************************************************/
-psF32 psModelEval(psModel *model, psImage *image, psS32 col, psS32 row)
-{
-    PS_PTR_CHECK_NULL(image, false);
-    PS_PTR_CHECK_NULL(model, false);
-    PS_PTR_CHECK_NULL(model->params, false);
+
+// XXX EAM : I have made this a public function
+// XXX EAM : this now uses a pmModel as the input
+// XXX EAM : it was using src->type to find the model, not model->type
+psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)
+{
+    PS_ASSERT_PTR_NON_NULL(image, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+    PS_ASSERT_PTR_NON_NULL(model->params, false);
 
     // Allocate the x coordinate structure and convert row/col to image space.
@@ -1449 +1329 @@
     x->data.F32[1] = (psF32) (row + image->row0);
     psF32 tmpF;
-    psModelFunc modelFunc;
-
-    modelFunc = psModelFunc_GetFunction (model->type);
+    pmModelFunc modelFunc;
+
+    modelFunc = pmModelFunc_GetFunction (model->type);
     tmpF = modelFunc (NULL, model->params, x);
     psFree(x);
@@ -1491 +1371 @@
     }
 
-    // XXX EAM : i changed this to match psModelEval above, but see
+    // XXX EAM : i changed this to match pmModelEval above, but see
     // XXX EAM   the note below in pmSourceContour
-    psF32 oldValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow);
+    psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
     if (oldValue == level) {
         return(((psF32) (subCol + source->pixels->col0)));
@@ -1514 +1394 @@
 
     while (subCol != lastColumn) {
-        psF32 newValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow);
+        psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
         if (oldValue == level) {
             return((psF32) (subCol + source->pixels->col0));
@@ -1674 +1554 @@
     //           tests below could be conditions (!NULL)
 
-    psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
 
     psVector *params = model->params;
@@ -1694 +1574 @@
     }
     if (count <  nParams + 1) {
-      psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
-      return(false);
+        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+        return(false);
     }
 
@@ -1722 +1602 @@
 
     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
-                            PM_SOURCE_FIT_MODEL_TOLERANCE);
+                                                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;
-      }
+        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;
+        }
     }       
 
@@ -1782 +1662 @@
 bool p_pmSourceAddOrSubModel(psImage *image,
                              psImage *mask,
-                             psModel *model,
+                             pmModel *model,
                              bool center,
                              psS32 flag)
 {
   
-  // XXX EAM : convert to ASSERTS
-  // PS_PTR_CHECK_NULL(model, false);
-  // PS_IMAGE_CHECK_NULL(image, false);
-  // PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+    PS_ASSERT_IMAGE_NON_NULL(image, false);
+    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
 
     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
     psVector *params = model->params;
-    psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     psS32 imageCol;
     psS32 imageRow;
@@ -1842 +1721 @@
 bool pmSourceAddModel(psImage *image,
                       psImage *mask,
-                      psModel *model,
+                      pmModel *model,
                       bool center)
 {
@@ -1852 +1731 @@
 bool pmSourceSubModel(psImage *image,
                       psImage *mask,
-                      psModel *model,
+                      pmModel *model,
                       bool center)
 {

trunk/psphot/src/pmObjects_EAM.h

@@ -1 @@
-
+/** @file  pmObjects.h
+ *
+ *  This file will ...
+ *
+ *  @author GLG, MHPCC
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-07 03:33:01 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+#if !defined(PM_OBJECTS_H)
+#define PM_OBJECTS_H
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include<stdio.h>
+#include<math.h>
+#include "pslib.h"
+
+// XXX EAM : this function should be added to psLib
+int          psImageCountPixelMask (psImage *mask, psU8 value);
+
+// XXX EAM : bug : missing from psLib/*.h
+bool psMinimizeGaussNewtonDelta (psVector *delta,
+                                 const psVector *params,
+                                 const psVector *paramMask,
+                                 const psArray  *x,
+                                 const psVector *y,
+                                 const psVector *yErr,
+                                 psMinimizeLMChi2Func func);
+
+// XXX EAM : switch to psPoly when ready
+psF64           Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y);
+
+// XXX EAM : psEllipse needs to be be included in psLib
+# include "psEllipse.h"
+
+/** pmPeakType
+ * 
+ *  A peak pixel may have several features which may be determined when the
+ *  peak is found or measured. These are specified by the pmPeakType enum.
+ *  PM_PEAK_LONE represents a single pixel which is higher than its 8 immediate
+ *  neighbors.  The PM_PEAK_EDGE represents a peak pixel which touching the image
+ *  edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific
+ *  number of neighbors at the same value, within some tolarence:
+ * 
+ */
+typedef enum {
+    PM_PEAK_LONE,                       ///< Isolated peak.
+    PM_PEAK_EDGE,                       ///< Peak on edge.
+    PM_PEAK_FLAT,                       ///< Peak has equal-value neighbors.
+    PM_PEAK_UNDEF                       ///< Undefined.
+} pmPeakType;
+
+/** pmPeak data structure
+ *  
+ *  
+ *  
+ */
+typedef struct
+{
+    int x;                              ///< X-coordinate of peak pixel.
+    int y;                              ///< Y-coordinate of peak pixel.
+    float counts;                       ///< Value of peak pixel (above sky?).
+    pmPeakType class;                   ///< Description of peak.
+}
+pmPeak;
+
+/** pmMoments data structure
+ *  
+ *  
+ *  
+ */
+typedef struct
+{
+    float x;                            ///< X-coord of centroid.
+    float y;                            ///< Y-coord of centroid.
+    float Sx;                           ///< x-second moment.
+    float Sy;                           ///< y-second moment.
+    float Sxy;                          ///< xy cross moment.
+    float Sum;                          ///< Pixel sum above sky (background).
+    float Peak;                         ///< Peak counts above sky.
+    float Sky;                          ///< Sky level (background).
+    float SN;                           ///< approx signal-to-noise
+    int nPixels;                        ///< Number of pixels used.
+} pmMoments;
+
+typedef int pmModelType;
+
+/** pmPSFClump data structure
+ *  
+ *  
+ *  
+ */
 typedef struct
 {
@@ -8 +106 @@
 } pmPSFClump;
 
-typedef struct
-{
-    float x;                            ///< X-coord of centroid.
-    float y;                            ///< Y-coord of centroid.
-    float Sx;                           ///< x-second moment.
-    float Sy;                           ///< y-second moment.
-    float Sxy;                          ///< xy cross moment.
-    float Sum;                          ///< Pixel sum above sky (background).
-    float Peak;                         ///< Peak counts above sky.
-    float Sky;                          ///< Sky level (background).
-    float SN;
-    int nPixels;                        ///< Number of pixels used.
-}
-pmMoments;
-
-typedef psS32 pmModelType;
-
-typedef struct
-{
-    pmModelType type;       ///< Model to be used.
-    psVector *params;       ///< Paramater values.
-    psVector *dparams;      ///< Parameter errors.
-    float radius;           ///< fit radius actually used
-    float chisq;            ///< Fit chi-squared.
-    int nDOF;             ///< number of degrees of freedom
-    int nIter;            ///< number of iterations to reach min
+/** pmModel data structure
+ *  
+ *  
+ *  
+ */
+typedef struct
+{
+    pmModelType type;                   ///< Model to be used.
+    psVector *params;                   ///< Paramater values.
+    psVector *dparams;                  ///< Parameter errors.
+    float chisq;                        ///< Fit chi-squared.
+    int nDOF;                           ///< number of degrees of freedom
+    int nIter;                          ///< number of iterations to reach min
+    float radius;                       ///< fit radius actually used
 }
 pmModel;
 
+/** pmSourceType enumeration
+ *  
+ *  
+ *  
+ */
 typedef enum {
     PS_SOURCE_DEFECT,
@@ -58 +149 @@
 } pmSourceType;
 
-typedef struct
-{
-    pmPeak *peak;           ///< Description of peak pixel.
-    psImage *pixels;        ///< Rectangular region including object pixels.
-    psImage *noise;         ///< Mask which marks pixels associated with objects.
-    psImage *mask;          ///< Mask which marks pixels associated with objects.
-    pmMoments *moments;     ///< Basic moments measure for the object.
-    pmModel *modelPSF;      ///< PSF Model fit (parameters and type)
-    pmModel *modelFLT;      ///< FLT Model fit (parameters and type).
-    pmSourceType type;      ///< Best identification of object.
+/** pmSource data structure
+ *  
+ *  This source has the capacity for several types of measurements. The
+ *  simplest measurement of a source is the location and flux of the peak pixel
+ *  associated with the source:
+ *  
+ */
+typedef struct
+{
+    pmPeak *peak;                       ///< Description of peak pixel.
+    psImage *pixels;                    ///< Rectangular region including object pixels.
+    psImage *noise;                     ///< Mask which marks pixels associated with objects.
+    psImage *mask;                      ///< Mask which marks pixels associated with objects.
+    pmMoments *moments;                 ///< Basic moments measure for the object.
+    pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
+    pmModel *modelFLT;                  ///< FLT Model fit (parameters and type).
+    pmSourceType type;                  ///< Best identification of object.
 }
 pmSource;
 
-// XXX EAM function to return pmModel
-typedef psMinimizeLMChi2Func pmModelFunc;
-typedef psF64 (*pmModelFlux)(const psVector *params);
-typedef bool  (*pmModelGuessFunc)(pmModel *model, pmSource *source);
-typedef bool  (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
-typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
-
-typedef bool  (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
-typedef bool  (*pmModelFitStatusFunc)(pmModel *model);
-
-// XXX EAM structure to carry model options
-typedef struct {
-  char *name;
-  int nParams;
-  pmModelFunc      modelFunc;
-  pmModelFlux      modelFlux;
-  pmModelRadius    modelRadius;
-  pmModelLimits        modelLimits;
-  pmModelGuessFunc modelGuessFunc;
-  pmModelFromPSFFunc modelFromPSFFunc;
-  pmModelFitStatusFunc modelFitStatusFunc;
-} pmModelGroup;
-
-// XXX EAM : strucures to define elliptical shape parameters
-typedef struct {
-  double major;
-  double minor;
-  double theta;
-} EllipseAxes;
-
-typedef struct {
-  double x2;
-  double y2;
-  double xy;
-} EllipseMoments;
-
-typedef struct {
-  double sx;
-  double sy;
-  double sxy;
-} EllipseShape;
-
-EllipseAxes EllipseMomentsToAxes (EllipseMoments moments);
-EllipseShape EllipseAxesToShape (EllipseAxes axes);
-EllipseAxes EllipseShapeToAxes (EllipseShape shape);
+/** pmPSF data structure
+ *  
+ *  
+ *  
+ */
+typedef struct
+{
+    pmModelType type;                   ///< PSF Model in use
+    psArray *params;                    ///< Model parameters (psPolynomial2D)
+    float chisq;                        ///< PSF goodness statistic
+    int nPSFstars;                      ///< number of stars used to measure PSF
+}
+pmPSF;
+
+pmPeak *pmPeakAlloc(
+    int x,                              ///< Row-coordinate in image space
+    int y,                              ///< Col-coordinate in image space
+    float counts,                       ///< The value of the peak pixel
+    pmPeakType class                    ///< The type of peak pixel
+);
+
+pmMoments *pmMomentsAlloc();
+pmModel *pmModelAlloc(pmModelType type);
+pmSource  *pmSourceAlloc();
+
+/******************************************************************************
+pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector
+above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
+the location (x value) of all peaks.
+*****************************************************************************/
+psVector *pmFindVectorPeaks(
+    const psVector *vector,             ///< The input vector (float)
+    float threshold                     ///< Threshold above which to find a peak
+);
+
+/******************************************************************************
+pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
+above the given threshold.  Returns a psList containing the location (x/y
+value) of all peaks.
+*****************************************************************************/
+psArray *pmFindImagePeaks(
+    const psImage *image,               ///< The input image where peaks will be found (float)
+    float threshold                     ///< Threshold above which to find a peak
+);
+
+/******************************************************************************
+pmCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have
+a peak value above the given maximum, or fall outside the valid region.
+*****************************************************************************/
+psList *pmCullPeaks(
+    psList *peaks,                      ///< The psList of peaks to be culled
+    float maxValue,                     ///< Cull peaks above this value
+    const psRegion *valid               ///< Cull peaks otside this psRegion
+);
+
+/******************************************************************************
+pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius):
+
+*****************************************************************************/
+pmSource *pmSourceLocalSky(
+    const psImage *image,               ///< The input image (float)
+    const pmPeak *peak,                 ///< The peak for which the pmSource struct is created.
+    psStatsOptions statsOptions,        ///< The statistic used in calculating the background sky
+    float innerRadius,                  ///< The inner radius of the suqare annulus for calculating sky
+    float outerRadius                   ///< The outer radius of the suqare annulus for calculating sky
+);
+
+/******************************************************************************
+ *****************************************************************************/
+bool pmSourceMoments(
+    pmSource *source,                   ///< The input pmSource for which moments will be computed
+    float radius                        ///< Use a circle of pixels around the peak
+);
 
 /******************************************************************************
 pmSourcePSFClump(pmArray *source, psMetaDeta *metadata): find the source PSF clump
- *****************************************************************************/
-pmPSFClump pmSourcePSFClump(psArray *source, ///< The input psSource
-                            psMetadata *metadata ///< Contains classification parameters
-    );
+*****************************************************************************/
+pmPSFClump pmSourcePSFClump(
+    psArray *source,                    ///< The input pmSource
+    psMetadata *metadata                ///< Contains classification parameters
+);
 
 /******************************************************************************
 pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the
 source classification.
- *****************************************************************************/
-bool pmSourceRoughClass(psArray *source, ///< The input psSource
-                            psMetadata *metadata, ///< Contains classification parameters
-                            pmPSFClump clump ///< Statistics about the PSF clump
-    );
-
-pmModel *pmSourceModelGuess(pmSource *source, ///< The input psSource
-                                    pmModelType model ///< The type of model to be created.
-    );
-
-// XXX EAM : added utility functions
-pmModelType        pmModelSetType (char *name);
-pmModelFunc        pmModelFunc_GetFunction (pmModelType type);
-pmModelFlux        pmModelFlux_GetFunction (pmModelType type);
-pmModelGuessFunc   pmModelGuessFunc_GetFunction (pmModelType type);
-pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type);
-pmModelRadius      pmModelRadius_GetFunction (pmModelType type);
-char                  *pmModelGetType (pmModelType type);
-psS32                  pmModelParameterCount (pmModelType type);
-pmModelLimits          pmModelLimits_GetFunction (pmModelType type);
-pmModelFitStatusFunc   pmModelFitStatusFunc_GetFunction (pmModelType type);
-int                    pmCountSatPixels (psImage *image);
-
-pmMoments *pmMomentsAlloc();
-pmModel *pmModelAlloc(pmModelType type);
-pmSource *pmSourceAlloc();
+*****************************************************************************/
+bool pmSourceRoughClass(
+    psArray *source,                    ///< The input pmSource
+    psMetadata *metadata,               ///< Contains classification parameters
+    pmPSFClump clump                    ///< Statistics about the PSF clump
+);
+
+/******************************************************************************
+pmSourceSetPixelCircle(source, image, radius)
+*****************************************************************************/
+bool pmSourceSetPixelsCircle(
+    pmSource *source,                   ///< The input pmSource
+    const psImage *image,               ///< The input image (float)
+    float radius                        ///< The radius of the circle
+);
+
+/******************************************************************************
+ *****************************************************************************/
+pmModel *pmSourceModelGuess(
+    pmSource *source,                   ///< The input pmSource
+    pmModelType model                   ///< The type of model to be created.
+);
+
+/******************************************************************************
+ *****************************************************************************/
+typedef enum {
+    PS_CONTOUR_CRUDE,
+} pmContourType;
+
+psArray *pmSourceContour(
+    pmSource *source,                   ///< The input pmSource
+    const psImage *image,               ///< The input image (float) (this arg should be removed)
+    float level,                        ///< The level of the contour
+    pmContourType mode                  ///< Currently this must be PS_CONTOUR_CRUDE
+);
+
+/******************************************************************************
+ *****************************************************************************/
+bool pmSourceFitModel(
+    pmSource *source,                   ///< The input pmSource
+    pmModel *model,                     ///< model to be fitted
+    const bool PSF                      ///< Treat model as PSF or FLT?
+);
+
+/******************************************************************************
+ *****************************************************************************/
+bool pmSourceAddModel(
+    psImage *image,                     ///< The output image (float)
+    psImage *mask,                      ///< The image pixel mask (valid == 0)
+    pmModel *model,                     ///< The input pmModel
+    bool center                         ///< A boolean flag that determines whether pixels are centered
+);
+
+/******************************************************************************
+ *****************************************************************************/
+bool pmSourceSubModel(
+    psImage *image,                     ///< The output image (float)
+    psImage *mask,                      ///< The image pixel mask (valid == 0)
+    pmModel *model,                     ///< The input pmModel
+    bool center                         ///< A boolean flag that determines whether pixels are centered
+);
+
+/**
+ * 
+ *  The object model functions are defined to allow for the flexible addition
+ *  of new object models. Every object model, with parameters represented by
+ *  pmModel, has an associated set of functions which provide necessary support
+ *  operations. A set of abstract functions allow the programmer to select the
+ *  approriate function or property for a specific named object model.
+ * 
+ */
+
+/**
+ * 
+ *  This function is the model chi-square minimization function for this model.
+ * 
+ */
+typedef psMinimizeLMChi2Func pmModelFunc;
+
+
+/**
+ * 
+ * This function returns the integrated flux for the given model parameters.
+ */
+typedef psF64 (*pmModelFlux)(const psVector *params);
+
+
+/**
+ * 
+ *  This function returns the radius at which the given model and parameters
+ *  achieves the given flux.
+ * 
+ */
+typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
+
+/**
+ * 
+ *  This function sets the model parameter limits vectors for the given model
+ * 
+ */
+typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
+
+/**
+ * 
+ *  This function provides the model guess parameters based on the details of
+ *   the given source.
+ * 
+ */
+typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
+
+
+/**
+ * 
+ *  This function constructs the PSF model for the given source based on the
+ *  supplied psf and the FLT model for the object.
+ * 
+ */
+typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
+
+
+/**
+ * 
+ *  This function returns the success / failure status of the given model fit
+ * 
+ */
+typedef bool (*pmModelFitStatusFunc)(pmModel *model);
+
+/**
+ * 
+ *  Each of the function types above has a corresponding function which returns
+ *  the function given the model type:
+ * 
+ */
+pmModelFunc          pmModelFunc_GetFunction (pmModelType type);
+pmModelFlux          pmModelFlux_GetFunction (pmModelType type);
+pmModelRadius        pmModelRadius_GetFunction (pmModelType type);
+pmModelLimits        pmModelLimits_GetFunction (pmModelType type);
+pmModelGuessFunc     pmModelGuessFunc_GetFunction (pmModelType type);
+pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction (pmModelType type);
+pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
+
+// pmModelGroup utility functions
+int                  pmModelParameterCount (pmModelType type);
+char                *pmModelGetType (pmModelType type);
+pmModelType          pmModelSetType (char *name);
+
+// structure to carry model group functions
+typedef struct {
+    char *name;
+    int nParams;
+    pmModelFunc          modelFunc;
+    pmModelFlux          modelFlux;
+    pmModelRadius        modelRadius;
+    pmModelLimits        modelLimits;
+    pmModelGuessFunc     modelGuessFunc;
+    pmModelFromPSFFunc   modelFromPSFFunc;
+    pmModelFitStatusFunc modelFitStatusFunc;
+} pmModelGroup;
+
+#endif

trunk/psphot/src/pmPeaksSigmaLimit.c

@@ -16 +16 @@
 
     psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32);
-    psImageSmooth (smooth, SIGMA, NSIGMA);
+    psImageSmooth_EAM (smooth, SIGMA, NSIGMA);
     psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
 

trunk/psphot/src/pmSourceUtils.c

@@ -7 +7 @@
                           psF32 Radius)
 {
-    psRegion *srcRegion;
+    psRegion srcRegion;
 
     // Grab a subimage of the original image of size (2 * outerRadius).
-    srcRegion = psRegionSquare (x, y, Radius);
-    srcRegion = psRegionForImage (srcRegion, imdata->image, srcRegion);
+    srcRegion = psRegionForSquare (x, y, Radius);
+    srcRegion = psRegionForImage (imdata->image, &srcRegion);
 
-    //  use these when psImageSubset is updated
-    //  psImage *subImage = psImageSubset(imdata->image, srcRegion);
-    //  psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
-    //  psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
+    mySource->pixels = psImageSubset(imdata->image, srcRegion);
+    mySource->noise  = psImageSubset(imdata->noise, srcRegion);
+    mySource->mask   = psImageSubset(imdata->mask,  srcRegion);
 
-    mySource->pixels = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
-    mySource->noise  = psImageSubset(imdata->noise, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
-    mySource->mask   = psImageSubset(imdata->mask,  srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
-
-    psFree (srcRegion);
     return(mySource);
 }

trunk/psphot/src/psImageData.c

@@ -20 +20 @@
     imdata->mask = mask;
 
-    p_psMemSetDeallocator(imdata, (psFreeFcn) psImageDataFree);
+    psMemSetDeallocator(imdata, (psFreeFunc) psImageDataFree);
     return (imdata);
 }

trunk/psphot/src/psLibUtils.c

@@ -1 @@
-# include "psphot.h"
+# include <strings.h>  // for strncasecmp
+# include <pslib.h>
+
+// XXX EAM : this is NOT included in the C99 headers ??
+FILE *fdopen(int fildes, const char *mode);
 
 // XXX EAM : these utility functions should be added back into PSLib
@@ -29 +33 @@
     if (timers == NULL) timers = psHashAlloc (16);
 
-    start = psTimeGetTime (PS_TIME_UTC);
+    start = psTimeGetNow (PS_TIME_UTC);
     psHashAdd (timers, name, start);
     psFree (start); 
@@ -47 +51 @@
     if (start == NULL) return (0);
 
-    mark = psTimeGetTime (PS_TIME_UTC);
+    mark = psTimeGetNow (PS_TIME_UTC);
     delta = psTimeDelta (mark, start);
     psFree (mark);
@@ -53 +57 @@
 
     return (delta);
+}
+
+// find the location of the specified argument
+int psArgumentGet (int argc, char **argv, char *arg) {
+
+    int i;
+
+    for (i = 0; i < argc; i++) {
+        if (!strcmp(argv[i], arg))
+            return (i);
+    }
+  
+    return ((int)NULL);
+}
+
+// remove the specified argument (by location)
+int psArgumentRemove (int N, int *argc, char **argv) {
+
+    int i;
+
+    if ((N != (int)NULL) && (N != 0)) {
+        (*argc)--;
+        for (i = N; i < *argc; i++) {
+            argv[i] = argv[i+1];
+        }
+    }
+
+    return (N);
 }
 
@@ -114 +146 @@
 }
 
-// find the location of the specified argument
-int psArgumentGet (int argc, char **argv, char *arg) {
-
-    int i;
-
-    for (i = 0; i < argc; i++) {
-        if (!strcmp(argv[i], arg))
-            return (i);
-    }
-  
-    return ((int)NULL);
-}
-
-// remove the specified argument (by location)
-int psArgumentRemove (int N, int *argc, char **argv) {
-
-    int i;
-
-    if ((N != (int)NULL) && (N != 0)) {
-        (*argc)--;
-        for (i = N; i < *argc; i++) {
-            argv[i] = argv[i+1];
-        }
-    }
-
-    return (N);
-}
-
 // alternate implementation of this function from pmObjects.c
+// XXX EAM : move this to pmObjects_EAM.c
 psVector *psGetRowVectorFromImage(psImage *image,
                                   psU32 row)
 {
-    PS_IMAGE_CHECK_NULL(image, NULL);
-    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
 
     psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
@@ -201 +206 @@
 }
 
-void psImageSmooth (psImage *image, float sigma, float Nsigma) {
+// XXX EAM : this is now in psLib
+void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma) {
 
     int Nx, Ny, Npixel, Nrange;
@@ -314 +320 @@
 }           
 
+// count number of pixels with given mask value
+int psImageCountPixelMask (psImage *mask, psU8 value) 
+{
+    int Npixels = 0;
+  
+    for (int i = 0; i < mask->numRows; i++) {
+        for (int j = 0; j < mask->numCols; j++) {
+            if (mask->data.U8[i][j] & value) {
+                Npixels ++;
+            }
+        }
+    }
+    return (Npixels);
+}
+
 // define a square region centered on the given coordinate
+// XXX EAM : this is now in psLib
+# if (0)
 psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) {
     psRegion *region;
@@ -321 +344 @@
     return (region);
 }
+# endif
 
 // set actual region based on image parameters:
@@ -327 +351 @@
 //     frame, which means the negative values should subtract from Nx,Ny of
 //     the parent, not the child.  but, we don't carry the dimensions of the 
-//     parent in the psImage container.  for now, us the child Nx,Ny
-// force range to be on this subimage 
-// XXX EAM : this needs to be changes to use psRegion rather than psRegion*
+//     parent in the psImage container.  for now, use the child Nx,Ny
+//     force range to be on this subimage 
+// XXX EAM : this needs to be changed to use psRegion rather than psRegion*
+// XXX EAM : this is now in psLib
+# if (0)
 psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) {
     
@@ -356 +382 @@
     return (out);
 }
+# endif
 
 // mask the area contained by the region
 // the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
 void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
 
@@ -375 +404 @@
     }
 }
+# endif
 
 // mask the area not contained by the region
 // the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
 void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
 
@@ -396 +428 @@
     }
 }
+# endif
 
 // mask the area contained by the region
 // the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
 void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
 
@@ -419 +454 @@
     }
 }
+# endif
 
 // mask the area contained by the region
 // the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
 void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
 
@@ -442 +480 @@
     }
 }
+# endif
 
 psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) {
@@ -455 +494 @@
   return (out);
 }
+
+// XXX EAM a utility function
+bool p_psVectorPrintRow (int fd, psVector *a, char *name)
+{
+
+    FILE *f;
+    f = fdopen(fd, "a+");
+    fprintf (f, "vector: %s\n", name);
+
+    for (int i = 0; i < a[0].n; i++) {
+        fprintf (f, "%f  ", p_psVectorGetElementF64(a, i));
+    }
+    fprintf (f, "\n");
+    fclose(f);
+    return (true);
+}
+// XXX EAM is the use of fdopen here a good way to do this?

trunk/psphot/src/psLine.c

@@ -14 +14 @@
   psLine *line;
   line = psAlloc (sizeof(psLine));
-  p_psMemSetDeallocator(line, (psFreeFcn) psLineFree);
+  psMemSetDeallocator(line, (psFreeFunc) psLineFree);
 
   line->Nline = 0;

trunk/psphot/src/psMinimize.c

@@ -13 +13 @@
                             psMinimizeLMChi2Func func)
 {
-    PS_PTR_CHECK_NULL(min, NULL);
-    PS_VECTOR_CHECK_NULL(params, NULL);
-    PS_VECTOR_CHECK_EMPTY(params, NULL);
-    PS_PTR_CHECK_NULL(x, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_EMPTY(y, NULL);
-    PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL);
-    PS_PTR_CHECK_NULL(func, NULL);
+    // PS_PTR_CHECK_NULL(min, NULL);
+    // PS_VECTOR_CHECK_NULL(params, NULL);
+    // PS_VECTOR_CHECK_EMPTY(params, NULL);
+    // PS_PTR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_NULL(y, NULL);
+    // PS_VECTOR_CHECK_EMPTY(y, NULL);
+    // PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL);
+    // PS_PTR_CHECK_NULL(func, NULL);
 
     // this function has test and current values for several things

trunk/psphot/src/psModulesUtils.c

@@ -11 +11 @@
     psImage *image = source->pixels;
     psImage *mask  = source->mask;
-    psPeak *peak  = source->peak;
-    psRegion *srcRegion;
-
-    srcRegion = psRegionSquare (peak->x, peak->y, Radius);
-    srcRegion = psRegionForImage (srcRegion, mask, srcRegion);
-    psImageMaskRegion (mask, srcRegion, OR, 0x80);
-
+    pmPeak *peak  = source->peak;
+    psRegion srcRegion;
+
+    // XXX EAM : psRegionXXX funcs need value not ptr
+
+    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, 0x7f);
+    psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED);
 
     psF64 tmpF64;
     p_psGetStatValue(myStats, &tmpF64);
-
     psFree (myStats);
-    psFree (srcRegion);
 
     if (isnan(tmpF64)) return (false);
@@ -40 +40 @@
                              const bool PSF)
 {
-    PS_PTR_CHECK_NULL(source, false);
-    PS_PTR_CHECK_NULL(source->moments, false);
-    PS_PTR_CHECK_NULL(source->peak, false);
-    PS_PTR_CHECK_NULL(source->pixels, false);
-    PS_PTR_CHECK_NULL(source->mask, false);
-    PS_PTR_CHECK_NULL(source->noise, false);
+    // PS_PTR_CHECK_NULL(source, false);
+    // PS_PTR_CHECK_NULL(source->moments, false);
+    // PS_PTR_CHECK_NULL(source->peak, false);
+    // PS_PTR_CHECK_NULL(source->pixels, false);
+    // PS_PTR_CHECK_NULL(source->mask, false);
+    // PS_PTR_CHECK_NULL(source->noise, false);
     psBool fitStatus = true;
     psBool onPic     = true;
@@ -188 +188 @@
  
 XXX: mask values?
+
+XXX EAM : this version clips input pixels on S/N
 *****************************************************************************/
 # define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
@@ -194 +196 @@
                          psF32 radius)
 {
-    PS_PTR_CHECK_NULL(source, NULL);
-    PS_PTR_CHECK_NULL(source->peak, NULL);
-    PS_PTR_CHECK_NULL(source->pixels, NULL);
-    PS_PTR_CHECK_NULL(source->noise, NULL);
+    // PS_PTR_CHECK_NULL(source, NULL);
+    // PS_PTR_CHECK_NULL(source->peak, NULL);
+    // PS_PTR_CHECK_NULL(source->pixels, NULL);
+    // PS_PTR_CHECK_NULL(source->noise, NULL);
     PS_FLOAT_COMPARE(0.0, radius, NULL);
 

trunk/psphot/src/psPolynomials.c

@@ -1 +1 @@
 # include "psphot.h"
+
+// XXX EAM : this file defines alternate versions of the polynomial fitting
+//           functions.  these differ from psLib v.0.5.0 in the definition 
+//           of nOrder (here we use nOrder in the mathematical sense, not
+//           like IRAF.  ie, a linear function is of nOrder 1)
+//           psLib v 0.7.0 should be in sync with out defintion of nOrder
+//           If so, we can drop these functions and use the psLib versions.
 
 // write out the terms of the given 1D polynomial
@@ -9 +16 @@
 }    
 
-psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly)
+psF64 Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly)
 {
     psS32 loop_x = 0;
-    psF32 polySum = 0.0;
-    psF32 xSum = 1.0;
+    psF64 polySum = 0.0;
+    psF64 xSum = 1.0;
 
     for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
@@ -28 +35 @@
                                        const psVector *x)
 {
-    PS_POLY_CHECK_NULL(myPoly, NULL);
-    PS_VECTOR_CHECK_NULL(x, NULL);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    // PS_POLY_CHECK_NULL(myPoly, NULL);
+    // PS_VECTOR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
 
     psVector *tmp;
@@ -56 +63 @@
                                   psPolynomialType type)
 {
-    PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL);
+    // PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL);
 
     psS32 i = 0;
@@ -65 +72 @@
     newPoly->type = type;
     newPoly->n = nOrder;
-    newPoly->coeff = (psF32 *)psAlloc(nTerm * sizeof(psF32));
-    newPoly->coeffErr = (psF32 *)psAlloc(nTerm * sizeof(psF32));
+    newPoly->coeff = (psF64 *)psAlloc(nTerm * sizeof(psF64));
+    newPoly->coeffErr = (psF64 *)psAlloc(nTerm * sizeof(psF64));
     newPoly->mask = (psU8 *)psAlloc(nTerm * sizeof(psU8));
     for (i = 0; i < nTerm; i++) {
@@ -73 +80 @@
         newPoly->mask[i] = 0;
     }
-    p_psMemSetDeallocator(newPoly, (psFreeFcn) psPolynomial1DFree);
+    psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial1DFree);
     return(newPoly);
 }
@@ -120 +127 @@
     // dump minutiae
 # ifndef PS_NO_TRACE
+# if (0)
+    // XXX EAM : change from FILE to fd breaks this code:
     if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) {
       FILE *f = psTraceGetDestination ();
@@ -127 +136 @@
       }
     }
+# endif
 # endif
       
@@ -214 +224 @@
 
 {
-    PS_POLY_CHECK_NULL(myPoly, NULL);
-    PS_VECTOR_CHECK_NULL(x, NULL);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+    // PS_POLY_CHECK_NULL(myPoly, NULL);
+    // PS_VECTOR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    // PS_VECTOR_CHECK_NULL(y, NULL);
+    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
 
     psVector *tmp;
@@ -246 +256 @@
 
 {
-    PS_POLY_CHECK_NULL(myPoly, NULL);
-    PS_VECTOR_CHECK_NULL(x, NULL);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+    // PS_POLY_CHECK_NULL(myPoly, NULL);
+    // PS_VECTOR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    // PS_VECTOR_CHECK_NULL(y, NULL);
+    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
 
     psVector *tmp;
@@ -291 +301 @@
                                     psPolynomialType type)
 {
-    PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
-    PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
+    // PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
+    // PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
 
     psS32 x = 0;
@@ -305 +315 @@
     newPoly->nY = nYorder;
 
-    newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
-    newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
-    newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
+    newPoly->coeff = (psF64 **)psAlloc(nXterm * sizeof(psF64 *));
+    newPoly->coeffErr = (psF64 **)psAlloc(nXterm * sizeof(psF64 *));
+    newPoly->mask = (char **)psAlloc(nXterm * sizeof(char *));
+    // XXX EAM : this is an error in the definition of the polynomial mask 
+    //           (should be psU8 not char)
     for (x = 0; x < nXterm; x++) {
-        newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
-        newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
-        newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
+        newPoly->coeff[x] = (psF64 *)psAlloc(nYterm * sizeof(psF64));
+        newPoly->coeffErr[x] = (psF64 *)psAlloc(nYterm * sizeof(psF64));
+        newPoly->mask[x] = (char *)psAlloc(nYterm * sizeof(char));
     }
     for (x = 0; x < nXterm; x++) {
@@ -320 +332 @@
         }
     }
-    p_psMemSetDeallocator(newPoly, (psFreeFcn) psPolynomial2DFree);
+    psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial2DFree);
     return(newPoly);
 }
@@ -510 +522 @@
                             const psVector* dz)
 {
-    PS_VECTOR_CHECK_NULL(mask, NULL);
-    PS_VECTOR_CHECK_NULL(x, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_NULL(z, NULL);
-    PS_VECTOR_CHECK_NULL(dz, NULL);
+    // PS_VECTOR_CHECK_NULL(mask, NULL);
+    // PS_VECTOR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_NULL(y, NULL);
+    // PS_VECTOR_CHECK_NULL(z, NULL);
+    // PS_VECTOR_CHECK_NULL(dz, NULL);
 
     psVector *zFit   = NULL;
@@ -548 +560 @@
 // XXX EAM : VectorFitPolynomial2DOrd and Polynomial2DEvalVector require different types (F32 vs F64)
 
-# if (0) // moved to psLib
-// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
-psF32 Polynomial2DEval(const psPolynomial2D* myPoly,
-                       psF32 x,
-                       psF32 y)
-{
-    PS_POLY_CHECK_NULL(myPoly, NAN);
+psF64 Polynomial2DEval(const psPolynomial2D* myPoly,
+                       psF64 x,
+                       psF64 y)
+{
+    // PS_POLY_CHECK_NULL(myPoly, NAN);
 
     psS32 loop_x = 0;
     psS32 loop_y = 0;
-    psF32 polySum = 0.0;
-    psF32 xSum = 1.0;
-    psF32 ySum = 1.0;
+    psF64 polySum = 0.0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
 
     // XXX EAM : nX is order, not nTerms
@@ -577 +587 @@
     return(polySum);
 }
-# endif
-

trunk/psphot/src/psphot.h

@@ +1 @@
+# include <stdio.h>
+# include <strings.h>  // for strcasecmp
+# include <unistd.h>   // for unlink
 # include <pslib.h>
-// # include <pmObjects.h>
-// my additions and modifications:
 # include "pmObjects_EAM.h"
-# include <string.h>
-# include <strings.h>
-# include <unistd.h>
-# include <stdlib.h>
-# include <math.h>
-
-// # define  M_PI 3.14159265358979323846264338328      /* pi */
 
 typedef struct {
-  psImage *image;
-  psImage *mask;
-  psImage *noise;
-  psMetadata *header;
+    psImage *image;
+    psImage *mask;
+    psImage *noise;
+    psMetadata *header;
 } psImageData;
 
 // data to test a given PSF model type
 typedef struct {
-  pmModelType modelType;
-  pmPSF      *psf;
-  psArray    *sources;      // pointers to the original sources
-  psArray    *modelFLT;     // model fits, floating parameters 
-  psArray    *modelPSF;     // model fits, PSF parameters
-  psVector   *mask;
-  psVector   *metric;
-  psVector   *fitMag;
-  float       ApResid;
-  float       dApResid;
-  float       skyBias;
+    pmModelType modelType;
+    pmPSF      *psf;
+    psArray    *sources;      // pointers to the original sources
+    psArray    *modelFLT;     // model fits, floating parameters 
+    psArray    *modelPSF;     // model fits, PSF parameters
+    psVector   *mask;
+    psVector   *metric;
+    psVector   *fitMag;
+    float       ApResid;
+    float       dApResid;
+    float       skyBias;
 } pmPSF_Test;
 
 // structure to carry a dynamic string
 typedef struct {
-  int NLINE;
-  int Nline;
-  char *line;
+    int NLINE;
+    int Nline;
+    char *line;
 } psLine;
 
@@ -43 +37 @@
     PSPHOT_MASK_INVALID   = 0x01,
     PSPHOT_MASK_SATURATED = 0x02,
-    PSPHOT_MASK_KEEP      = 0x08,
+    PSPHOT_MASK_MARKED    = 0x08,
 } MaskValues;
-
-// used by mask operations
-# define AND true
-# define OR false
 
 # define psMemCopy(A)(psMemIncrRefCounter((A)))
@@ -116 +106 @@
 psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
 psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type);
-//void       psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue);
-//void       psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue);
-//void       psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue);
-//void       psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue);
 psVector    *psGetRowVectorFromImage(psImage *image, psU32 row);
+int          psImageCountPixelMask (psImage *mask, psU8 value);
 
 // basic image functions
 bool         psImageInit (psImage *image,...);
-//void       psImageSmooth (psImage *image, float sigma, float Nsigma);
-//psRegion    *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
-psRegion    *psRegionSquare (psF32 x, psF32 y, psF32 radius);
+void         psImageSmooth_EAM (psImage *image, float sigma, float Nsigma);
 int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
 
@@ -135 +120 @@
 
 // polynomial functions
-psF32           Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
+psF64           Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y);
 psImage        *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
 psPolynomial2D *VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
@@ -143 +128 @@
 psPolynomial2D *RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
 psVector       *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
+
 psVector       *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
 void            psPolynomial1DDump (psPolynomial1D *poly);
-psF32           Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly);
+psF64           Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly);
 psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
 psPolynomial1D *Polynomial1DAlloc(psS32 nOrder, psPolynomialType type);
 psPolynomial1D *VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
 
-// fitsource
-bool            onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y);
-psMetadata     *fs_args (int *argc, char **argv);
-int             fs_usage ();
-bool            DumpImage (psImage *image, char *filename);
+// XXX EAM : bug : missing from psLib/*.h
+bool psMinimizeGaussNewtonDelta (psVector *delta,
+                                 const psVector *params,
+                                 const psVector *paramMask,
+                                 const psArray  *x,
+                                 const psVector *y,
+                                 const psVector *yErr,
+                                 psMinimizeLMChi2Func func);
+
+// not included in the .h file?
+bool p_psVectorPrint (
+    int fd,                            ///< output file descriptor
+    psVector *a,                       ///< vector to print
+    char *name                         ///< name of vector (for title)
+);
+
+bool p_psVectorPrintRow (int fd, psVector *a, char *name);

trunk/psphot/src/psphotApplyPSF.c

@@ -66 +66 @@
 
         // fit PSF model (set/unset the pixel mask)
-        psImageKeepCircle (source->mask, x, y, model->radius, "or", PSPHOT_MASK_KEEP);
+        psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
         status = pmSourceFitModel (source, model, true);
-        psImageKeepCircle (source->mask, x, y, model->radius, "and", ~PSPHOT_MASK_KEEP);
+        psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
         if (!status || (model->params->data.F32[1] < 0)) {
           psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);

trunk/psphot/src/psphotArguments.c

@@ -43 +43 @@
   psMetadata *config = psMetadataAlloc ();
   psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_META_MULTI, "folder for psf model entries", NULL);
-  config = psMetadataParseConfig (config, &Nfail, argv[3], FALSE);
+  config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
   fprintf (stderr, "loaded config...\n");
 

trunk/psphot/src/psphotFitGalaxies.c

@@ -50 +50 @@
 
         // recalculate the source moments using the larger galaxy moments radius
-        status = pmSourceMoments (source, GAL_MOMENTS_RAD);
+        // XXX EAM : 0.5.0 code used _EAM version
+        // status = pmSourceMoments (source, GAL_MOMENTS_RAD);
+        status = pmSourceMoments_EAM (source, GAL_MOMENTS_RAD);
         if (!status) {
           source->type = PS_SOURCE_DROP_GALAXY;  // better choice?
@@ -74 +76 @@
 
         // fit FLT (not PSF) model (set/unset the pixel mask)
-        psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP);
+        psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
+        // XXX EAM : 0.5.0 code used _EAM version
+        // status = pmSourceFitModel (source, model, false);
         status = pmSourceFitModel (source, model, false);
-        psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP);
+        psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
         if (!status) {
           // if the fit fails, we need to change the classification

trunk/psphot/src/psphotOutput.c

@@ -372 +372 @@
 
 // write the peaks to an output file
-bool psPeaksWriteText (psArray *sources, char *filename) {
+bool pmPeaksWriteText (psArray *sources, char *filename) {
 
     int i;
     FILE *f;
-    psPeak *peak;
+    pmPeak *peak;
         
     f = fopen (filename, "w");
     if (f == NULL) {
-        psLogMsg ("psPeaksWriteText", 3, "can't open output file for peaks%s\n", filename);
-        return false;
-    }
-
-    for (i = 0; i < sources->n; i++) {
-        peak = (psPeak *) sources->data[i];
+        psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename);
+        return false;
+    }
+
+    for (i = 0; i < sources->n; i++) {
+        peak = (pmPeak *) sources->data[i];
         if (peak == NULL) continue;
         fprintf (f, "%5d %5d  %7.1f\n", 
@@ -638 +638 @@
     unlink (filename);
     psFits *fits = psFitsAlloc (filename);
-    psFitsWriteImage (fits, header, image, 0, NULL);
+    psFitsWriteImage (fits, header, image, 0);
     psFree (fits);
     return (TRUE);

trunk/psphot/src/psphotSetup.c

@@ -75 +75 @@
     psRegion keep = psRegionSet (XMIN, XMAX, YMIN, YMAX);
     keep = psRegionForImage (image, &keep);
-    psImageKeepRegion (mask, &keep, OR, PSPHOT_MASK_INVALID);
+    psImageKeepRegion (mask, &keep, "OR", PSPHOT_MASK_INVALID);
     // XXX EAM : psRegionForImage should not take psRegion *keep
 

trunk/psphot/src/psphotSourceStats.c

@@ -35 +35 @@
 
         // measure basic source moments
-        status = pmSourceMoments (source, RADIUS);
+        // XXX EAM : choose between these two versions
+        status = pmSourceMoments_EAM (source, RADIUS);
         if (!status) {
           psFree (source);

trunk/psphot/src/pspsf.c

@@ -31 +31 @@
     }
 
-    p_psMemSetDeallocator(psf, (psFreeFcn) pmPSFFree);
+    psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
     return(psf);
 }
@@ -74 +74 @@
     }   
 
-    p_psMemSetDeallocator(test, (psFreeFcn) pmPSF_TestFree);
+    psMemSetDeallocator(test, (psFreeFunc) pmPSF_TestFree);
     return (test);
 }
@@ -106 +106 @@
         // set temporary object mask and fit object
         // fit model as FLT, not PSF
-        psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
+        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
         status = pmSourceFitModel (source, model, false);
-        psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
 
         // exclude the poor fits
@@ -142 +142 @@
         y = source->peak->y;
 
-        psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
+        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
         status = pmSourceFitModel (source, modelPSF, true);
 
@@ -168 +168 @@
 
     next_source:
-        psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
 
     }