Index: trunk/psModules/src/objects/pmPeaks.c
===================================================================
--- trunk/psModules/src/objects/pmPeaks.c	(revision 9881)
+++ trunk/psModules/src/objects/pmPeaks.c	(revision 10048)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-11-07 09:07:42 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-11-17 23:00:31 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -24,28 +24,49 @@
 
 /******************************************************************************
-myListAddPeak(): A private function which allocates a psArray, if the list
-argument is NULL, otherwise it adds the peak to that list.
-XXX EAM : changed the output to psArray
-XXX EAM : Switched row, col args
-XXX EAM : NOTE: this was changed in the call, so the new code is consistent
-*****************************************************************************/
-static psArray *myListAddPeak(psArray *list,
-                              psS32 row,
-                              psS32 col,
-                              psF32 counts,
-                              pmPeakType type)
-{
-    psTrace("psModules.objects", 4, "---- %s() begin ----\n", __func__);
-    pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
-
-    if (list == NULL) {
-        list = psArrayAllocEmpty(100);
-    }
-    psArrayAdd(list, 100, tmpPeak);
-    psFree (tmpPeak);
-    // XXX EAM : is this free appropriate?  (does psArrayAdd increment memory counter?)
-
-    psTrace("psModules.objects", 4, "---- %s() end ----\n", __func__);
-    return(list);
+AddPeak(): A private function which allocates a psArray, if the peaks
+argument is NULL, otherwise it adds the peak to that array.
+XXX EAM : row,col now refer to image coords, NOT parent (since this is private) 
+XXX EAM : now also calculates fractional peak positions from 3x3 bicube region
+*****************************************************************************/
+static psArray *AddPeak(psArray *peaks,
+                        const psImage *image,
+                        psS32 row,
+                        psS32 col,
+                        pmPeakType type)
+{
+    psTrace(__func__, 5, "---- begin ----\n");
+
+    if (peaks == NULL) {
+        peaks = psArrayAllocEmpty(100);
+    }
+
+    // the peak position is in parent coordinates
+    pmPeak *peak = pmPeakAlloc(col + image->col0, row + image->row0, image->data.F32[row][col], type);
+
+    // measure fractional peak position using the 3x3 bicube fit
+
+    // ix,iy must land on image with 1 pixel border
+    int ix = PS_MAX (PS_MIN (col, image->numCols - 2), 1);
+    int iy = PS_MAX (PS_MIN (row, image->numRows - 2), 1);
+
+    // calculate peak position relative to ix,iy
+    psPolynomial2D *bicube = psImageBicubeFit (image, ix + image->col0, iy + image->row0);
+    psPlane min = psImageBicubeMin (bicube);
+    psFree (bicube);
+
+    // if min point is too deviant, use the peak value
+    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
+        peak->xf = min.x + ix + image->col0;
+        peak->yf = min.y + iy + image->row0;
+    } else {
+        peak->xf = ix;
+        peak->yf = iy;
+    }
+
+    psArrayAdd(peaks, 100, peak);
+    psFree (peak);
+
+    psTrace(__func__, 5, "---- end ----\n");
+    return(peaks);
 }
 
@@ -306,10 +327,5 @@
     psArray *list = NULL;
 
-    psU32 col0 = image->col0;
-    psU32 row0 = image->row0;
-
-    //
     // Find peaks in row 0 only.
-    //
     row = 0;
     tmpRow = getRowVectorFromImage((psImage *) image, row);
@@ -319,7 +335,5 @@
     for (psU32 i = 0 ; i < row1->n ; i++ ) {
         col = row1->data.U32[i];
-        //
-        // Determine if pixel (0,0) is a peak.
-        //
+        // is pixel (0,0) is a peak?
         if (col == 0) {
             if ( (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
@@ -328,5 +342,5 @@
 
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = AddPeak(list, image, row, col, PM_PEAK_EDGE);
                 }
             }
@@ -338,5 +352,5 @@
                     (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = AddPeak(list, image, row, col, PM_PEAK_EDGE);
                 }
             }
@@ -347,5 +361,5 @@
                     (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = AddPeak(list, image, row, col, PM_PEAK_EDGE);
                 }
             }
@@ -386,5 +400,5 @@
                         (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
                     myType = PM_PEAK_EDGE;
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
+                    list = AddPeak(list, image, row, col, myType);
                 }
             } else if (col < (image->numCols - 1)) {
@@ -421,13 +435,5 @@
                         }
 
-                        list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
-
-                        # if (0)
-
-                            psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
-                        psPlane min = psImageBicubeMin (bicube);
-                        peak->xf = min.x;
-                        peak->yf = min.y;
-                        # endif
+                        list = AddPeak(list, image, row, col, myType);
 
                     }
@@ -442,5 +448,5 @@
                         (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
                     myType = PM_PEAK_EDGE;
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
+                    list = AddPeak(list, image, row, col, myType);
                 }
             } else {
@@ -466,5 +472,5 @@
                     (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = AddPeak(list, image, row, col, PM_PEAK_EDGE);
                 }
             }
@@ -476,5 +482,5 @@
                     (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = AddPeak(list, image, row, col, PM_PEAK_EDGE);
                 }
             }
@@ -485,5 +491,5 @@
                     (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
                 if (image->data.F32[row][col] > threshold) {
-                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
+                    list = AddPeak(list, image, row, col, PM_PEAK_EDGE);
                 }
             }